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I.  INTRODUCTION* 


A.  GENERAL  BACKGROUND 

Space  is  an  extremely  hostile  environment.  Whether  in  free  space,  orbiting  the 
earth,  or  on  the  surface  of  foreign  worlds,  mankind  is  ill  equipped  to  work  or  even 
survive  for  long  periods  of  time.  Yet  mankind  is  driven  by  its  natural  curiosity  to 
explore  this  vast  frontier.  Not  only  is  mankind  attempting  to  explore  space  but 
plans  are  underway  to  build  habitats  in  space  and  to  live  off  earth  in  the  near 
future.  There  are  bold  concepts  being  developed  for  factories  in  space,  huge 
colonies  of  people  orbiting  between  the  Moon  and  Earth,  and  manned  landings  on 
Mars  within  the  next  fifty  years.  However,  the  weak  link  in  all  these  ambitious 
plans  is  mankind  itself.  Tied  to  bulky  space  suits  and  physiologically  limited  by 
the  amount  of  exposure  to  microgravity  and  radiation  that  can  be  withstood, 
mankind  is  barely  capable  of  doing  the  large  scale  manual  labor  required  to  bring 
these  dreams  to  fruition.  Consequently,  there  are  those  that  feel  that  if  mankind  is 
to  ever  truly  become  a  spacefaring  being  it  will  have  to  be  done  through  a 
symbiotic  relationship  with  intelligent  machines  or  robots. 

In  this  relationship,  robots  will  be  utilized  to  perform  those  functions  for 
which  it  is  not  economically  feasible  to  task  humans  with.  These  functions  will 
include  manual  labor  such  as  the  construction  and  maintenance  of  large  scale 
space  structures  where  transportation  costs  will  make  it  prohibitive  to  employ 


people  as  common  metal  workers.  Robots  will  also  be  delegated  the  tasks  too 
boring  or  dangerous  for  mankind  to  perform  such  as  doing  long  term  exploration 
missions  on  the  surface  of  the  planets  which  mankind  may  wish  to  investigate. 
Certainly,  mankind  will  soon  return  to  the  surface  of  the  Moon  and  will  someday 
land  on  Mars  but  it  will  be  robots  which  will  be  tasked  to  spend  the  boring  weeks 
and  months  exploring  the  surface. 

For  its  part,  mankind  will  perform  those  tasks  which  robots  are  ill  equipped 
to  do.  These  include  analyzing  the  data  and  samples  returned  by  robotic 
explorers,  doing  precision  manual  labor  such  as  maintaining  electrical  equipment, 
and  other  tasks  requiring  a  high  degree  of  cognitive  or  reasoning  capability. 
Mankind  will  be  the  thinking  part  of  the  team  and  the  robots  will  be  the  muscle. 

But  just  as  mankind  will,  on  occasion,  be  required  to  use  muscle  in  the 
conquest  of  space,  so  will  robots  be  required  to  exhibit  intelligent  behavior.  An 
autonomous  vehicle  exploring  the  surface  of  Mars  would  be  an  economic  disaster, 
considering  transportation  and  fuel  costs,  if  every  time  the  autonomous  vehicle 
came  to  a  cliff  it  would  just  walk  right  over  the  edge  or  if  the  vehicle  ran  into 
every  hazard  in  its  path  instead  of  taking  the  required  action  to  avoid  that 
hazard.  However,  today  the  technology  required  to  produce  the  necessary 
intelligent  behavior  in  robots  is  just  beginning  to  be  developed.  These  technologies 
include  artificial  intelligence  programming,  walking  technology,  and  vision  sensors 
to  name  but  a  few. 
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The  purpose  of  this  thesis  will  be  to  investigate,  by  simulation,  a  new 
algorithm  for  computer  vision  for  autonomous  vehicle  navigation  and  path 
planning.  This  algorithm  will  attempt  to  classify  terrain  by  type  rather  than 
specifically  identify  individual  terrain  items. 

B.  ORGANIZATION 

A  discussion  of  the  recent  history  in  computer  vision  research  is  presented  in 
Chapter  2.  This  chapter  includes  presentations  on  some  of  the  basic  work  in  vision 
research,  a  number  of  the  different  sensor  types  used  to  implement  vision,  and  the 
most  recent  w-ork  utilizing  state  of  the  art  testbeds.  Finally,  the  background  for 
the  terrain  classification  algorithm  implemented  in  this  simulation  is  addressed. 

In  Chapter  3.  a  detailed  mathematical  development  of  the  formulas  used  in 
the  simulation  is  given.  A  complete  understanding  of  the  mathematics  is  required 
to  understand  the  limitations  on  the  algorithm  and  some  of  the  requirements 
necessary  in  implementing  it  in  hardware. 

Chapter  4  provides  a  description  of  the  program  written  for  the  simulation. 
This  includes  a  brief  discussion  of  the  LISP  progamming  language  and  bitmap 
graphics,  both  of  which  are  important  to  this  simulation.  Finally,  the 
methodology  for  displaying  the  symbolic  information  generated  by  the  program  is 
discussed. 

The  symbolic  results  from  the  simulation  are  presented  in  Chapter  5.  The 
simulation  produces  images  of  the  terrain  scenes  analyzed  and  symbolic  maps  of 


each  scene  as  generated  by  the  algorithm-  Each  scene  is  analyzed  under  a  number 
of  different  conditions.  These  conditions  include  both  noiseless  as  well  as  noisy 
data  and  varying  internal  thresholds.  The  effect  of  each  of  these  variables  is 
briefly  discussed. 

Finally,  Chapter  6  summarizes  the  contributions  made  by  this  research  and 
discusses  some  of  the  possible  areas  of  further  research  based  on  this  work.  Three 
appendices,  a  secondary  mathematical  development  of  the  formulas  as  based  on 
matrix  algebra,  the  graphics  routines  utilized,  and  the  simulation  as  it  is  written 
in  LISP,  are  included  after  Chapter  6. 
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II.  SURVEY  OF  PREVIOUS  WORK 


A.  INTRODUCTION 

To  make  a  land  vehicle  truly  autonomous,  the  principles  of  robotics  have  to 
be  applied  in  its  design.  These  principles  include  giving  the  vehicle  the  ability  to 
perceptive  and  understand  the  local  environment  from  sensor  data  and  the 
artificial  intelligence  to  plan  it’s  actions  as  a  result  of  this  understanding.  To  one 
degree  or  another,  all  animal  life  shares  a  common  set  of  environmental  sensors 
and  so,  intuitively,  those  sensors  are  the  ones  that  roboticists  attempt  to 
incorporate  in  their  designs. 

The  set  of  natural  senses  includes  touch,  taste,  hearing,  smell,  and  sight. 
Taste  and  smell  are  primarily  chemically-based  sensors.  Touch  and  hearing  are 
pressure  related  phenomenon.  Finally,  the  ability  to  sense  temperature  remotely 
and  to  see  are  essentially  passive  electromagnetic  capabilities  that  operate  in  two 
different  wavelength  regimes.  Although  all  the  senses  must  be  fully  operational  for 
a  person  to  feel  comfortable  in  their  ability  to  completely  sense  the  local 
environment,  the  one  natural  sense  that  seems  to  be  the  primary  sensor,  at  least 
for  human  beings,  is  vision. 

While  all  the  senses  are  normally  required,  a  person  is  able  to  adequately 
perform  basic  functions  even  with  the  loss  of  most  of  the  other  senses  so  long  as 
the  vision  system  is  not  impaired.  Accordingly,  the  ability  to  simulate  vision  is 


one  of  the  essential  requirements  for  a  high  performance  autonomous  vehicle  [Ref. 
l:pp.  79-121].  This  ability  to  simulate  vision  has  been,  over  the  last  few  decades, 
one  of  the  central  areas  of  robotics  research. 

Until  recently,  most  vision  research  has  had  to  contend  with  extracting 
information  from  two-dimensional  data  only.  This  is  because  the  conventional 
mechanism  for  passively  receiving  electromagnetic  energy  has  been  a  camera  of 
some  sort.  Whether  these  cameras  have  been  standard  TV  cameras,  CCD  devices, 
or  sensitive  IR  receivers,  their  output  is  two-dimensional  in  nature.  To  most 
researchers,  this  is  acceptable  since  that  is  the  way  humans  see.  However,  it  now 
appears  that  the  human  vision  mechanism  is  extremely  sophisticated  in  the  way 
that  it  processes  this  two-dimensional  data  to  extract  the  third  dimension  of  range 
from  a  scene.  The  massive  parallelism  used  by  the  brain  in  its  signal  processing 
algorithms  will  be  very  difficult  to  duplicate,  at  least  in  the  near  future  [Ref.  2:pp. 
180-190].  Clearly,  though,  range  is  an  important  type  of  information  for  any 
autonomous  vehicle.  Consequently,  work  in  robotic  image  processing  has 
progressed  from  two-dimensional  image  segmentation  and  pattern  recognition  into 
active  three-dimensional  image  sensing.  [Ref.  3:pp.  206-220  ;  Ref.  4:pp.  3-12  ;  Ref. 
5:pp.  1-15] 

The  purpose  of  this  chapter  is  to  introduce  some  of  the  specifics  of  range 
sensing,  the  techniques  of  scene  description,  and  finally,  to  examine  previous 
approaches  to  feature  extraction  based  on  three-dimensional  imaging. 


B.  APPROACHES  TO  ACQUIRING  RANGE  IMAGES 
1.  Conventional  Radar 

The  earliest  technology  for  directly  acquiring  a  range  image  was 
conventional  radar.  The  term  RADAR  is  an  acronym  for  RAdio  Detection  And 
Ranging.  As  the  acronym  implies,  radar  emits  radio  frequency  energy  and  uses  the 
information  contained  in  the  reflected  beam  to  derive  information  regarding  the 
target.  The  beam  of  radio  frequency  energy  can  be  in  pulse  form  or  continuous 
wave  form  [Ref.  6:pp.  1-40,359-391]. 

For  conventional  radar,  the  usual  format  is  the  pulse  form  in  which  the 
radar  beam  is  pulsed  on  and  off  at  a  specified  rate.  This  rate  is  called  the  pulse 
repetition  frequency  ,  PRF.  The  length  of  time  that  the  pulse  is  on  is  called  the 
pulse  width  ,  PW.  These  two  parameters,  along  with  the  operating  frequency  and 
aperture  dimensions  of  the  transmitting/receiving  antenna,  define  most  of  the 
important  characteristics  of  the  radar. 

a.  Range  Measurement 

The  range  information  in  a  radar  system  is  derived  from  time  of  flight 
calculations.  The  pulsed  radar  beam  travels  at  the  speed  of  light,  C,  3  x  108 
meters  per  second.  By  measuring  the  time  from  when  the  beam  is  transmitted  to 
when  the  reflected  beam  returns,  and  dividing  that  time  by  two  to  account  for  the 
two  way  trip,  the  range  can  be  calculated  as: 


(2.1) 


Mt) 

2 


>:  C. 


If  this  was  all  there  was  to  measuring  range,  then  the  range  to  any  target  could  be 
determined  exactly. 

Unfortunately,  this  is  not  the  case.  Receiving  mechanisms  have  a 
limited  threshold  for  detecting  the  reflected  signal.  A  minimum  amount  of  energy 
has  be  be  received  before  the  detector  can  determine  that  the  energy  is  a  signal 
and  not  just  noise.  The  amount  of  time  after  a  receiver  starts  receiving  a  reflected 
signal  until  it  can  classify  the  signal  as  more  than  noise  is  a  function  of  the  power 
in  the  returned  signal,  since  energy  is  defined  as  power  times  time.  The  time  lost 
in  signal  determination  introduces  a  error  in  the  range  calculation.  As  the  range 
to  the  target  increases,  the  power  in  the  received  signal  is  reduced  as  a  function  of 
the  fourth  power  of  the  range,  so  the  time  from  the  instant  when  the  signal  is 
initially  received  till  when  it  can  be  classified  as  true  signal  increases  as  a  function 
of  range.  This  variation  in  power  with  distance  contributes  another  significant 
error  in  the  range  determination.  Finally,  few  radars  are  capable  of  actually 
determining  the  range  from  the  leading  edge  of  the  pulse  only.  Instead, 
information  is  derived  from  the  entire  length  of  the  pulse.  Therefore,  the  pulse 
width  is  commonly  used  as  a  lower  limit  for  the  range  resolution  of  the  radar 
beam.  For  example,  a  PW  of  .3  microseconds  has  an  inherent  range  resolution  of 
90  meters.  This  distance  is  also  considered  as  the  lower  limit  on  the  minimum 
range  that  a  radar  can  measure. 
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An  obvious  solution  to  the  problem  of  range  resolution  would  be  to 
reduce  the  PW  to  an  absolute  minimum.  To  a  degree,  this  is  in  fact  done. 
However,  minimizing  PW  is  not  without  its  own  set  of  limiting  problems.  The 
initial  power  of  the  transmitted  pulse  has  to  be  large  enough  so  that  the  energy  in 
the  returned  pulse  will  exceed  the  thresholding  limits  of  the  receiver.  So,  as  pulse 
width  is  reduced  to  increase  range  resolution,  the  power  of  the  RF  pulse  generator 
has  to  go  up  proportionately. 

The  PRF  also  has  a  significant  effect  on  range  measurement  with 
conventional  radar.  Generally,  the  antenna  that  is  used  to  transmit  the  high 
power  radar  beam  is  the  same  one  that  is  used  to  receive  the  extremely  weak 
reflected  signal.  Clearly,  the  radar  cannot  transmit  and  receive  at  the  same  time. 
If  there  is  even  a  small  impedance  mismatch  at  the  antenna/waveguide  interface, 
then  a  portion  of  the  radar  signal  is  reflected  back  to  the  transmitter/receiver.  If 
the  receiver  were  on  while  the  beam  was  being  transmitted,  the  energy  reflected 
by  the  interface  would,  at  a  minimum,  swamp  out  the  returning  signal  and.  in  the 
worst  case,  it  could  physically  destroy  the  extremely  sensitive  receiver. 

The  maximum  range  the  radar  can  detect  is  defined  by  the  PRF.  The 
on  and  off  time  of  the  transmitter  defines  one  cycle.  The  PRF  is  the  number  of 
on-off  cycles  per  second.  The  on  time  during  one  cycle  is  defined  by  the  PW,  so 
that: 
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off -time  = 


1 


-  PW  . 


(2.2) 


PRF 

During  the  off  time,  the  beam  has  to  travel  to  the  target  and  back  again  to  the 
receiver.  The  maximum  detectable  range  is.  then,  one  half  the  pulse  off-time  times 
the  speed  of  light. 

The  radar  parameters  that  directly  effect  the  ability  to  measure  range 
are,  in  summary,  pulse  width,  which  determines  power  and  range  resolution,  and 
pulse  repetition,  rate  which  determines  maximum  range, 
b.  Lateral  Measurement 

Raw  radar  spatial  measurements  are  usually  recorded  in  spherical 
coordinates.  Spherical  coordinates  are  defined  by  two  angles,  usually  theta  (  0  ) 
and  phi  (  4>  ),  and  range,  r.  Theta  is  an  azimuth  angle  measured  either  between  -ff 
and  ?r  or  0  and  2tt  radians  from  a  known  reference.  Phi  is  an  elevation  angle 
measured  between  -7r/2  and  7r/2  from  a  reference  elevation  (usually  the  horizon). 
Range  and  the  ability  to  resolve  range  has  been  discussed  in  the  previous 
paragraphs.  Measuring  elevation  and  azimuth  are  relatively  simple  in  radar.  The 
elevation  and  azimuth  are  simply  the  measured  position  of  the  antenna  at  the 
time  of  transmit /receive.  However,  resolving  differences  between  small  angular 
displacements  places  a  limit  on  the  ability  to  measure  physical  dimensions  and 
angular  resolution  is  a  function  of  the  beam  width. 

From  optical  diffraction  theory,  a  generalized  formula  for  determining 
the  beam  width,  and  therefore  the  angular  limit  of  resolution,  of  an 
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AS  1  22  *  D  (2.3) 

where  a  is  the  waveirngtt.  trie  radar  and  D  is  the  diameter  of  the  antenna  (or 
focusing  lens,  This  result  i*  m  radians  To  get  an  answer  in  terms  of  resolvable 
distance.  AS  is  multiplied  bv  the  range  to  the  target  i.e.. 

d  -  AS  Range  (2.4) 

As  an  example,  the  wavelength  for  a  1  ghz  radar  beam  is  10  cm. 
Assuming  a  1  meter  antenna  and  a  range  to  the  target  of  1000  meters,  the 
minimum  resolvable  separation  will  be  122  meters.  This  is  acceptable  for  target 
detection  and  a  rough  determination  of  azimuth,  but  for  an  autonomous  vehicle, 
this  is  not  very  useful. 

c.  Continuous  Wave  Radar 

Continuous  wave  radar  is  not  discussed  as  an  imaging  technology 
since  it  has  limited  capability  as  a  imaging  system.  It  is  generally  used  as  a  target 
illumination  system  or  as  a  doppler  moving  target  detector.  Range  detection  is 
limited  to  continuous  wave  radars  which  have  been  modulated  by  a  periodic  wave 
and  the  range  accuracy  of  such  systems  is  limited  [Ref.  6:.pp  359-390]. 


d.  Conclusions 


Radar  was  one  of  the  earliest  technologies  to  use  some  of  the  concepts 
of  range  imaging.  The  limitations  of  radar,  though,  which  are  inherent  in  the 
required  pulse  widths,  pulse  repetition  frequencies,  power  limitations,  and  beam 
width  due  to  frequency,  make  it  an  unattractive  solution  to  the  range  imaging 
requirements  of  a  terrestrial  autonomous  vehicle. 

2.  Infrared.  Millimeter  Wave  and  SAR 

Higher  frequency  (shorter  wavelength)  active  systems  and  advanced 
technolog}’  systems  offer  increased  accuracy  in  pointing  and  raiige  measurement. 
These  systems  include  infrared,  millimeter  wave  technology  and  synthetic 
aperture  radars.  At  present,  there  is  significant  research  into  this  kind  of  advanced 
technology.  However,  for  now,  their  capabilities  are  limited  in  scope  when  used  to 
analyze  natural  scenes.  Also,  the  equipment  is  generally  either  large,  fragile, 
unreliable,  or  expensive.  Accordingly,  the  military  is  the  main  user  of  these 
advanced  technologies  and  their  adaptation  to  civilian  uses  (such  as  robotics)  will 
take  time.  [Ref.  l:pp.  107] 

3.  Optical  Stereo 

Another  region  of  the  electromagnetic  spectrum  in  which  the  receiver 
technology  is  rather  well  developed,  though  less  sensitive,  is  the  visible  region.  A 
common  way  to  generate  range  information  in  this  regime  is  by  the  use  of  a 
technique  called  disparity  stereo  or  binocular  vision.  As  shown  in  Figure  1,  this 
technique  requires  two  optical  sensors  in  the  same  plane,  viewing  the  same  scene 


Binoc 


with  their  central  axes  parallel.  Range,  z.  to  the  target.  P.  is  measured  in  terms  of 


the  angles  .  02  and  the  baseline,  b.  between  the  two  receivers.  From  the  law  of 
sines 


Rx  R2 


b 

sin  (180° — 0j-02) 


(2.5) 


but 


The  theory  of  stereo  disparity  is  simple  and  straightforward.  However, 
implementing  it  is  a  difficult  procedure  and  subject  to  significant  error.  There  are 
two  major  techniques  for  implementing  a  depth  measurement  scheme  via 
binocular  vision.  They  are  by  passive  and  active  means, 
a.  Passive  Stereo  Ranging 

,  This  approach  is  very  similar  to  human  vision.  The  distance  to  a 

point  is  inversely  proportional  to  the  disparity  between  the  two  image  planes.  The 
ranging  becomes  more  accurate  the  farther  apart  the  two  cameras  are  or  the  closer 
the  target  is.  If  the  two  cameras  are  too  far  apart,  though,  the  requirement  for 
overlapping  fields  of  view'  may  become  a  difficult  prerequisite  to  meet.  Bringing 
the  two  cameras  too  close  together  makes  the  disparity  between  the  two  image 
planes  very  small  and  introduces  error  into  the  calculations. 

The  most  significant  problem  with  this  technique  is  that  a  given 
point.  P.  has  to  be  accurately  identified  in  both  images.  If  the  location  of  P  is 
displaced  slightly  due  to  thresholding  or  calibration  discrepancies  between  the  two 
cameras,  the  range  determination  will  be  in  error.  Moreover,  just  identifying  P  in 
both  image  planes  is  a  problem.  Since  camera  output  is  nothing  more  than  a  grey 
scale  code  of  each  pixel,  and  the  position  of  P  in  each  image  plane  will  not  be  the 
same,  all  of  the  pixels  in  each  image  have  to  be  correlated  to  each  pixel  in  the 
other  image.  Segmentation  of  the  image  planes  then  has  to  be  done  and  each 
segment  has  to  be  examined  for  areas  of  high  correlation.  Due  to  ambient  lighting 
variations  between  the  image  planes,  these  correlation  functions  may  be  slightly  in 


error  and  therefore  the  precise  location  of  known  points  in  each  image  that  are 
also  in  the  other  image  cannot  be  accurately  determined  [Ref.  ll:pp.  15-23]. 

This  technique  is  computationally  intensive  and  results  can  only  be 
considered  accurate  if  a  large  number  of  correlated  points  in  each  image  can  be 
determined  so  that  statistical  methods  of  error  reduction  can  be  applied.  However, 
its  advantages  lie  in  the  fact  that  it  is  most  similar  to  natural  ranging  methods 
and  it  is  not  power  intensive.  Illumination  of  the  area  of  interest  is  supplied  by 
nature,  not  by  the  autonomous  vehicle.  This  natural  illumination  is  also,  though, 
the  source  of  the  most  significant  error  since  it  cannot  be  preconditioned  (an 
example  would  be  imaging  of  a  scene  during  night  time  with  a  lack  of  significant 
ambient  light)  and  therefore  is  not  considered  well-behaved  [Ref.  7:pp.  150-151  ; 
Ref.  8:pp.  111-114]. 

b.  Active  Stereo  Ranging 

A  solution  to  the  stereo  image  point  correspondence  problem  is  to 

provide  an  easy  to  identify  point.  With  an  optical  sensor,  that  could  be  a  point  of 

light  significantly  brighter  than  the  ambient  light  background.  In  this  approach,  a 

spot  of  light  or  a  laser  beam  is  scanned  over  a  scene  and  the  spot  is  tracked.  The 

requisite  angular  information  is  measured  in  the  normal  way  and  range  is 

computed  similarly  [Ref.  12:pp.  134-146].  The  limiting  factors  in  this  ranging 

system  are  the  inherent  ability  to  measure  small  angles,  the  spot  size  of  the  beam, 

and  the  on  board  power  requirements  for  the  light  source.  Another  method  of 

active  optical  ranging  was  studied  for  a  number  of  years  at  Rensselaer  Polytechnic 
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Institute  [Ref.  13] .  This  method  did  not  depend  on  stereo  sensors  for  it's 
triangulation  to  the  points  of  interest  but  rather  an  array  of  directional 
photodiodes.  These  photodiodes  would  detect  a  laser  beam  reflection  from  the 
surrounding  terrain  which  had  been  illuminated  by  a  laser  offset  on  a  mast.  By  an 
appropriate  set  of  transformations,  the  elevation  and  range  to  the  point  reflecting 
the  laser  could  then  be  calculated. 

A  further  enhancement  is  to  project  an  entire  line  rather  than  a  single 
spot.  The  line  is  scanned  over  the  target  and  the  reflected  image  is  recorded  by  a 
camera.  Points  of  interest  can  then  be  designated  and  triangulation  methods  can 
be  applied  to  determine  ranges  in  the  scene.  Surfaces  that  are  nearly  parallel  to 
the  illuminating  plane,  though,  can  return  poor  results.  This  can  be  overcome  by 
projecting  two  mutually  orthogonal  series  of  scanning  planes  of  light.  Even 
further,  a  grid  can  be  projected  on  the  scene  so  that  scanning  hardware  is 
eliminated  [Ref.  14:pp.  150-151  :  Ref.  15:pp.  112  :  Ref.  10:pp  167-169], 

As  the  system  is  made  more  complicated,  the  processing  requirements 
increase  proportionately.  Also,  in  any  case,  the  inherent  accuracy  of  the  system  is 
dependent  on  the  ability  to  measure  small  angles  precisely,  just  as  before.  As  the 
range  to  the  points  of  interest  increases,  the  angles  decrease  and  accurate  range 
determination  becomes  more  difficult.  While  the  points  of  interest  are  more  easily 
determined,  this  advantage  is  not  without  it's  cost.  The  power  required  to 
illuminate  the  target  is  significant  and,  when  operating  in  an  environment  where 


humans  will  be  present,  the  distraction  of  illumination  by  an  artificial  grid  may 
not  be  acceptable  [Ref.  14:pp.  150-151  ;  Ref.  12:pp.  112  ;  Ref.  10:pp  167-169]. 

4.  Optical  Time  of  Flight 

Optical  time  of  flight  combines  many  of  the  best  characteristics  of  radar 
and  of  active  stereo  disparity  and  few  of  the  disadvantages.  Instead  of  using  radio 
frequency  electromagnetic  energy,  with  it’s  associated  long  wavelength  and  pulse 
width,  optical  time  of  flight  uses  visible  wavelength  or  near  IR  energy.  Because  of 
it’s  similarity  to  radar,  optical  time  of  flight  range  imaging  is  sometimes  called 
LIDAR.  This  acronym  stands  for  Light  Ranging  and  Detection. 

Though,  lidar  operates  in  or  near  the  visible  band  of  electromagnetic 
radiation  ,  it  can  be  characterized  by  the  same  parameters  as  radar.  It’s  beam 
width  is  calculated  exactly  as  is  a  radar  beam.  Using  Equation  (2.2)  and  the 
typical  operating  wavelength  of  a  GaAs  diode  laser  of  850  nanometers,  a  1 
millimeter  diameter  focusing  lens  will  produce  a  beam  width  of  1.04  milliradians 
[Ref.  16:pp.  156-167],  At  1000  meters  (well  beyond  typical  ranging  distances  for 
lidar)  this  beam  has  a  beam  width  of  1.0  meters,  much  smaller  than  the 
corresponding  radar  beam.  Attendant  with  the  narrow  beam  width,  the  lidar  has 
a  significantly  increased  angular  resolution  when  it  is  imaging  a  target  and  so  it 
can  generate  a  much  more  precise  set  of  data. 

Noise  in  the  target  area  in  the  form  of  ambient  light  is  a  major  source  of 
error  in  binocular  vision  systems,  both  active  and  passive,  as  is  radio  noise  a 
source  of  ranging  error  for  radar.  The  narrow  band-width  of  a  laser  beam  when 


compared  to  the  relatively  wide  band  in  which  the  power  of  natural  visible  light  is 
distributed,  reduces  the  noise  problem  enough  so  that,  when  measuring  distance 
by  time  of  flight,  very  narrow  pulse  widths  can  be  used.  For  example,  a  one 
milliwatt  helium-neon  laser  operating  at  633  nanometers  with  a  beam  width  of  1 
milliradian  and  a  band  width  of  0.2  nanometers  has  a  spectral  brightness  of 
1  x  10  lumens/ cm  -steradian-nm.  The  sun's  spectral  brightness  is  it’s  luminance 
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(1.5  x  10  lumens/ cm  -steradian)  divided  by  it's  band  width  (approximately  300 
nanometers).  This  is  about  500  lumens/cm  -steradian-nm  [Ref.  16:pp.  20-23].  The 
low  power  laser  is  over  5  orders  of  magnitude  brighter  than  the  sun,  the  major 
source  of  optical  noise.  Realization  of  this  advantage  is  possible  to  some  degree  by 
using  optical  filters  to  eliminate  all  reflected  light  except  that  light  in  the  lasers 
band  width  although  available  filter  bandwidths  are  many  orders  of  magnitude 
wider  than  laser  bandwidths  [Ref.  I7:pp.  55-59]. 

Actually  measuring  range  by  a  lidar  can  be  accomplished  by  either 
pulsing  the  laser  and  measuring  the  time  of  reflection,  as  per  radar,  or  by 
modulating  a  continuous  wave  laser  beam  and  measuring  the  phase  difference  of 
the  reflected  beam.  This  second  method  has  the  added  advantage  of  allowing 
doppler  processing  of  the  frequency  of  the  reflected  energy  to  derive  velocity 
information  on  the  target  [Ref.  10:pp.  170  ;  Ref.  14:pp.  151-152  :  Ref.  15:pp.  113- 
114]. 

Although  the  pulsed  systems  are  more  similar  to  radar,  adequate  range 
resolution  and  the  minimum  detectable  ranges  required  by  autonomous  vehicles 


necessitate  very  short  pulses.  Pulses  in  the  range  of  10‘s  of  picoseconds  will  give 
range  resolutions  in  the  millimeter  range  'Ref.  14:pp.  151-152  :  Ref.  15:pp.  113- 
114].  Picosecond  pulses  are  possible  using  Q  switching  techniques,  but  that 
requires  sophisticated  and  expensive  equipment  and  the  energy  in  such  a  small 
pulse  is  very  low.  Phase  modulated  CW  ranging  systems  offer  good  range 
resolution  and  minimum  range  with  less  sophisticated  and  sensitive  equipment. 
Modulation  frequencies  of  10  to  50  megahertz  have  been  used  to  measure  ranges 
with  an  accuracy  of  few  millimeters  over  a  range  of  less  than  a  meter  out  to  more 
than  10  meters  [Ref.  10:pp.  170  ;  Ref.  14:pp.  151-152  ;  Ref.  15:pp.  113-114]. 

Optical  time  of  flight  range  imaging  systems  are  a  relatively  new 
technology.  They  promise  good  range  resolution,  very  narrow  beam-widths,  good 
noise  performance,  and  reasonable  ranging  distances  for  the  local  navigation 
problem.  Unlike  stereo  optical  techniques,  the  range  resolution  is  not  a  function  of 
the  range  from  the  sensor.  Their  beam-width  is  superior  to  conventional  radar 
resulting  in  superior  lateral  resolution.  As  the  technology  matures,  they  may 
ultimately  prove  to  be  the  dominant  form  of  local  range  imaging  for  robotics. 

C.  SYMBOLIC  DESCRIPTION  OF  TERRAIN  FOR  AUTONOMOUS 
VEHICLES 

Once  the  data  from  an  imaging  system  has  been  collected,  the  next  step  for 
an  autonomous  vehicle  is  to  glean  information  from  it  effectively.  An  array  of  grey 
scale  pixels  or  x,  y.  z  values  is,  in  and  of  itself,  rather  useless.  The  generation  of 


knowledge  about  the  local  environment  from  images  is  a  two  step  process.  These 
two  steps  have  been  appropriately  titled  as  early  and  late  visual  processing  [Ref. 
9:pp.  87-160].  In  early  processing,  the  goal  is  to  generate  characteristic 
information  about  the  scene  that  the  late  visual  processing  modules  can  use.  Early 
processing  has  been  the  center  of  most  vision  research  conducted  until  recently. 
The  information  generated  by  early  processing  modules  is  a  line  drawing  termed 
the  intrinsic  image  [Ref.  9:pp.  95].  The  intrinsic  image  is  defined  as  the  grouping 
of  the  image  into  regions  corresponding  to  physical  surfaces.  This  segmentation  is 
often  done  by  identifying  edge  pixels  in  a  scene  (through  edge  identification 
algorithms  which  measure  brightness  gradients),  connecting  these  edge  pixels  with 
lines  (line  generation  algorithms)  and  then  identifying  the  physical  dimensions  of 
the  regions  contained  within  the  line  boundaries.  While  this  is  the  most  common 
form  of  segmentation  for  the  intrinsic  image,  there  are  as  many  different  ways  to 
segment  a  scene  as  there  are  people  conducting  research  in  machine  vision.  Each 
segmentation  scheme  amounts  to  a  new  concept  in  the  effort  to  discover  the 
inherent  (or  intrinsic)  knowledge  contained  in  the  raw  data  obtained  by  optical 
sensors. 

Once  the  intrinsic  image  has  been  defined,  by  whatever  scheme,  the  next  level 
in  visual  processing  is  called  late  image  processing  or  sometimes  image 
understanding  [Ref.  9:pp.  87-160].  It  is  with  the  late  image  processing  that 
artificial  intelligence  is  most  often  associated.  The  difference  between  simply 
building  and  connecting  lines  and  the  act  of  discerning  knowledge  from  those 


collections  of  lines  is  intelligence.  The  present  level  of  intelligence  programming  on 
computers  is  at  the  level  of  symbolic  computing.  An  intrinsic  image  is  reduced  to 
a  symbolic  description  (length,  width,  color,  height,  radius  of  curvature,  etc...) 
and  then  an  attempt  to  match  that  symbolic  description  with  a  predefined  set  of 
known  symbols  is  made.  If  a  good  match  is  made  with  symbol  X,  then  the 
intrinsic  image  can  be  assumed  to  be  another  X  or  very  similar  to  an  X. 
Obviously,  an  autonomous  vehicle  cannot  contain  enough  memory  to  store  a 
complete  set  of  symbols  for  every  possible  object  in  the  world,  and  if  it  could,  the 
processing  time  to  search  the  entire  set  of  symbols  would  be  prohibitive.  For  an 
autonomous  vehicle  to  navigate  about  a  local  environment,  though,  the  required 
symbolic  terrain  data  base  can  be  made  more  reasonable. 

1.  Object  Identification 

There  are  two  methods  of  image  processing  in  the  autonomous  vehicle 
local  navigation  problem.  These  methods  can  be  termed  as  object  identification 
and  region  classification  ,  respectively.  The  most  common  method  is  to  attempt 
to  identify  individual  objects  in  the  local  environment  first.  Then  the  path 
planner  routines  make  appropriate  decisions  based  on  this  knowledge  [Ref.  18:pp. 
223-235].  This  object  identification  can  be  either  on  an  extremely  sophisticated 
level  or  on  the  simplest  of  levels. 

A  very  simple  system  might  only  identify  objects  as  obstacle  or 
traversable  space.  This  is  the  most  elementary  level  at  which  vision  works.  The 
simplest  way  to  think  of  the  symbolic  description  of  an  obstacle  is  in  terms  of 


voxel's.  A  voxel  (VOlume  piXEL)  is  an  incremental  space  in  range  and  lateral 
displacement  (dx.  dy  in  cartesian  coordinates)  corresponding  to  the  minimum 
range  increment  and  horizontal  separation  detectable  by  the  vision  system.  This  is 
equivalent  to  a  "range  bin"  in  spherical  coordinates.  An  obstacle  will  be  then,  a 
voxel  with  enough  vertical  development  to  impede  the  motion  of  the  autonomous 
vehicle.  This  minimum  development  may  be  anything  greater  than  1  inch  for  a 
small  wheeled  vehicle,  for  instance.  (A  more  sophisticated  system  might  have  to 
concern  itself  with  vertical  development  down  from  the  ceiling  for  a  tall  vehicle. 
An  example  of  this  would  be  a  doorway.)  Voxels  that  are  not  obstacles  are  then 
classified  as  traversable.  Their  symbolic  description  would  be  characterized  as  a 
lack  of  assigned  attributes  [Ref.  18:pp.  223-235]. 

Obstacle  elements  that  are  adjacent  can  be  joined  to  each  other  to  create 
the  intrinsic  image  [Ref.  19:pp.  119-123].  Any  traversable  spaces  that  are 
completely  enclosed  by  obstacles  should  then  be  classified  as  obstacles  also.  For 
this  simple  vision  system,  after  the  scene  has  been  segmented  into  obstacle  and 
traversable  space,  the  autonomous  vehicle  need  go  no  further  in  object 
recognition.  The  symbolic  description  of  the  scene  would  be  groups  of  obstacles, 
each  one  specified  by  a  position  and  shape.  In  this  case,  the  object  identification 
routines  are  not  concerned  with  questions  about  the  physical  nature  of  the 
obstacles.  It  makes  no  difference  if  the  obstacle  is  a  desk  or  a  wall.  Indeed,  the 
obstacle  may  be  a  desk  against  a  wall,  for  all  the  system  cares,  and  it  will  be 


treated  the  same. 


For  most  autonomous  systems,  the  above  level  of  sophistication  is  not 
sufficient.  It  only  provides  a  vehicle  with  the  ability  to  avoid  collisions  with 
obstacles.  It  cannot  provide  the  vehicle  with  sufficient  knowledge  of  a  scene  for 
the  vehicle  to  identify  specific  objects  in  that  scene.  The  more  the  autonomous 
vehicle  is  required  to  exhibit  intelligent  behavior  in  it's  interaction  with  the  local 
environment,  the  more  the  symbolic  description  of  the  scene  will  have  to  be 
detailed.  Absolute  precision  in  the  symbolic  description  of  the  terrain  will  seldom, 
if  ever,  however,  be  obtainable.  This  is  due  to  the  problem  of  self  occlusion.  The 
reverse  side  of  an  object  is  hidden  from  imaging  systems  and  objects  in  front  of 
others  will  hide  detail.  This  has  lead  to  range  images  being  classified  as  2  1/2 
dimensional  versus  3  dimensional  [Ref.  10:pp  158  ;  Ref.  20:pp.  69].  Thus,  any 
symbolic  description  of  an  object  will  be  partial  at  best  and  will  cause  errors  when 
matching  it  with  known  object  symbolic  descriptions. 

Once  an  intrinsic  image  has  been  created,  inherent  characteristics  of  an 
object  may  be  derived  [Ref.  21:pp.  3-4].  Those  inherent  characteristics  comprise 
what  will  be  called  the  symbolic  description  of  that  object  [Ref.  9:pp.  143-160]. 
These  characteristics  will  consist  of  such  things  as  the  x,  y.  and  z  dimensions  of 
the  object,  for  instance.  Derivation  of  these  characteristics  will  allow  an  object's 
parameters,  such  as  its  volume  and  its  major  and  minor  axis,  to  be  calculated. 
Radius  of  curvature  is  also  a  characteristic  which  is  easily  derivable  from  the 
intrinsic  image.  The  choice  of  the  characteristic  a  system  attempts  to  derive  from 
the  intrinsic  image  is  driven  by  the  type  of  symbols  that  are  chosen  to  describe 


the  objects.  Some  algorithms  attempt  to  describe  objects  in  terms  general 
parallelograms  (x.‘  v.  and  z  coordinates)  while  others  use  generalized  cylinders 
(radius  of  curvature)  [Ref.  22:pp  629-637  :  Ref.  19:pp.  119-123].  These  two 
characteristics  have  been  chosen  only  for  the  purposes  of  illustration.  Algorithms 
have  been  utilized  which  measure  the  smoothness  of  a  surface,  moments  of  the 
shapes,  and  many  many  more  [Ref.  5:pp.  1-15]. 

The  derivation  of  all  of  the  above  types  of  information  would  be 
extremely  expensive  in  terms  of  computational  time.  For  this  reason,  most  vision 
research  attempts  to  identify  a  few  important  inherent  characteristics  to  analyze 
and  then  to  optimize  the  algorithms  used  to  produce  those  characteristics  or 
implement  the  routines  in  hardware.  If  the  number  of  individual  objects  a  system 
has  to  identify  is  small,  then  the  number  of  symbolic  descriptors  required  can  be 
kept  proportionately  small  and  the  computational  time  required  to  generate  them 
can  be  kept  manageable.  However,  for  real  world  applications,  most  environments 
that  an  autonomous  system  will  have  to  navigate  in  will  be  complicated  scenes. 
Consider  an  office  space,  a  factory  or  a  typical  home,  for  instance.  The  number  of 
individual  objects  present  in  the  environment  at  any  one  time  will  be  large  and. 
even  for  an  algorithm  that  generates  the  minimum  number  of  symbolic 
descriptors,  the  computational  time  required  to  fully  describe  all  the  objects  in  the 
scene  will  be  very  large. 
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2.  Region  Classification 


For  an  autonomous  vehicle  to  navigate  in  a  local  environment,  it  may  not 
be  necessary  for  the  system  to  actually  identify  each  element  in  the  scene.  For 
example,  in  an  exterior  environment,  it  is  not  necessary  for  a  robot  to  differentiate 
between  a  telephone  pole  and  a  tree.  It  is  enough  to  recognize  the  presence  of  an 
obstacle.  However,  the  simple  obstacle  recognition  systems  are  not  sufficient 
either.  An  autonomous  vehicle  has  the  ability  to  navigate  over  a  variety  of  terrain 
types  and  the  vehicle  may  have  to  make  a  decision  on  the  basis  of  a  cost  function 
of  which  type  of  terrain  to  travel  over.  An  example  would  be  a  scene  with  a  small 
hill  directly  in  the  vehicle's  path  and  flat  terrain  off  to  the  side.  The  vehicle  would 
have  to  decide  on  traveling  over  the  hill  or  altering  it’s  path  enough  to  follow  the 
flat  terrain.  F or  this  type  of  vision  another  approach  to  object  recognition,  termed 
region  classification,  has  recently  been  proposed  [Ref.  23:p.  1146). 

The  concept  of  region  classification  is  introduced  to  avoid  the  time 
consuming  and  complex  process  of  object  identification  by  using  topographic 
features  as  the  basis  for  an  autonomous  vehicle’s  navigation  decisions. 
Topographic  features  can  be  discerned  rather  rapidly  based  on  the  gradient  and 
second  directional  derivatives  around  an  individual  voxel.  These  topographic 
features  include  holes  and  poles,  valleys  and  ridges,  flat  and  inclined  (safe  and 
unsafe  inclined)  terrain,  and  saddles  [Ref.  24:pp.  22-24].  For  an  autonomous 
vehicle,  this  vision  technique  would  classify  both  a  tree  and  a  telephone  pole  as. 
simply,  a  pole.  This  eliminates  the  extra  time  required  to  distinguish,  on  the  basis 


of  more  complex  symbolic. descriptors,  between  the  two  objects.  Once  a  voxel  has 
been  classified,  a  symbolic  map  of  the  local  terrain  can  be  constructed  for  use  by 
the  path  planning  module.  This  process  will  be  more  completely  described  later  in 
this  thesis. 

D.  SPECIAL  REQUIREMENTS  OF  WALKING  VEHICLES 

An  autonomous  vehicle  may  be  described  by  a  number  of  different 
classification  parameters,  all  of  which  have  a  significant  effect  on  the  vehicle's 
performance.  The  mobility  of  a  vehicle  is  typically  classified  into  two  major  areas. 
These  are  wheeled  and  legged  systems  [Ref.  25:pp.  1-15, pp.  73-74].  The  wheeled 
systems  may  be  further  divided  into  conventional  wheeled  and  tracked  systems. 
Each  type  of  system  has  it’s  strong  points  and  it's  deficiencies. 

Wheeled  and  tracked  systems  have  a  preferred  direction  of  motion.  In  a  well 
defined  environment,  such  as  an  interior  space,  directional  preference  may  prove 
to  be  only  a  small  deficiency  or  even  an  advantage.  Having  a  preferred  direction 
limits  the  vision  requirements  of  the  system  because  the  system  can  safely 
concentrate  it’s  vision  in  the  forward  direction.  With  the  ability  to  retain  a  model 
of  the  space  it  passes  over,  when  reverse  motion  is  required,  due  to  an  error  in 
navigation,  a  reverse  vision  capability  is  not  an  absolute  requirement.  This 
reduction  in  the  vision  system  requirements  limits  the  complexity  of  the  system. 
Also,  by  virtue  of  the  fact  that  wheeled  systems  have  been  utilized  for  a  very  long 
time,  such  systems  are  simpler  to  design  and  build. 


However,  in  a  poorly  defined  environment,  such  as  natural  terrain,  a  wheeled 
vehicle  is  at  a  definite  disadvantage  when  compared  to  a  legged  vehicle  [Ref.  12:p. 
134].  Legged  vehicles  have  no  preferred  direction  of  motion.  They  may  translate 
sideways  with  as  much  ease  as  maneuvering  forward  or  backward.  This  allows  a 
significant  increase  in  maneuverability  of  the  vehicle.  Additionally,  whereas  a 
wheeled  vehicle  is  limited  in  its  step  crossing  ability  by  the  diameter  of  its  wheels 
and  the  bottoming  of  its  body,  a  legged  vehicle  is  only  limited  by  the  working 
envelope  of  its  legs.  This  working  envelope  is  defined  by  the  physical  length  of  the 
legs  and  the  requirement  to  maintain  the  center  of  gravity  within  the  support 
base  of  the  legs  in  physical  contact  with  the  ground  [Ref.  25:pp.  73-74] .  Clearly,  a 
properly  designed  legged  vehicle  will  have  a  much  greater  step  crossing  ability 
than  a  similar  wheeled  vehicle. 

This  increase  in  maneuverability  is  not  without  penalties,  though.  First, 
because  the  mechanics  of  designing  and  building  legged  systems  is  not  well 
understood,  systems  development  takes  extensive  time  and  requires  significant 
resources.  It  is  only  within  the  last  twenty  or  so  years  that  experiments  with  such 
systems  have  been  possible.  For  the  most  part,  these  vehicles  have  been  slow  and 
have  had  a  rough  ride.  Within  the  last  decade,  the  first  few  legged  systems 
starting  to  exhibit  an  ability  to  simulate  biological  motion  have  been  built  and 
demonstrated  [Ref  .18:pp.  7-15].  Research  in  this  area  has  been  slow,  but  recently 
has  been  showing  significant  progress. 
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Secondly,  the  inherent  increase  in  mobility  associated  with  legged  vehicles 
leads  to  another  penalty  by  putting  significant  demands  on  the  associated  vision 
systems.  Since  such  systems  can  move  sideways  as  well  as  forwards  and 
backwards,  a  vision  system  designed  to  work  with  an  independent  autonomous 
legged  system  should  be  able  to  preview  the  terrain  to  both  sides  as  well  as  fore 
and  aft.  Additionally,  depending  on  the  type  of  terrain  a  vehicle  is  crossing,  the 
vision  system  must  be  able  to  classify  surface  variations  of  much  smaller 
magnitude  than  presently  required.  Specifically,  with  the  level  of  sensing 
demonstrated  to  date  for  walking  vehicles,  computer-controlled  limb  cycling  in 
response  to  human  commanded  vehicle  direction  and  speed  is  sufficient  to 
maneuver  a  walking  vehicle  over  smooth  terrain  only.  However,  because  a  hole  as 
small  as  the  foot  pad  could  prove  to  be  a  potentially  destructive  obstacle,  a  vision 
system  must  be  able  to  detect  such  an  object  and  to  accurately  sense  its  position 
[Ref.  12:pp.  134-135).  Also,  in  order  to  traverse  the  rough  terrain  these  legged 
vehicles  are  envisioned  for.  they  must  be  able  to  "see"  individual  foot  holds  for 
the  navigation  and  control  system  to  utilize;  otherwise  their  mobility  will  only  be 
marginally  better  than  wheeled  systems. 

E.  APPROACHES  TO  FEATURE  EXTRACTION 
1.  Robot  Rover 

One  of  the  earliest  attempts  to  combine  vision  with  an  autonomous 
vehicle  was  the  Robot  Rover  project  by  Hans  P.  Moravec  at  the  Stanford  AI  lab 


[Ref.  11;.  This  vehicle  used  stereo  vision  to  generate  a  three  dimensional  view  of 
the  local  environment.  One  of  the  unique  aspects  of  this  stereo  vision  system  was 
that  it  was  implemented  with  a  single  camera.  The  camera  was  mounted  on  a  52 
centimeter  slider  bar  above  a  simple  mechanized  cart.  At  each  sampling  point  the 
cart  would  stop  and  the  camera  would  slide  from  left  to  right  taking  nine  pictures 
at  6.5  centimeter  intervals.  An  interest  operator  would  then  identify 
approximately  thirty  points  in  the  fifth  (middle)  image  for  a  correlator  module  to 
attempt  to  match  in  the  other  images.  Once  the  points  were  matched,  their  range 
could  be  calculated,  off  cart,  on  a  DEC  KL10  computer  conected  by  means  of  a 
radio  link,  using  the  apparent  parallax  between  images.  The  reason  for  using  nine 
images  was  to  attempt  to  reduce  the  effect  of  noise  in  generating  range  errors 
since  the  correlator  would  often  identify  features  incorrectly.  Specifically,  in 
Moravec’s  work,  the  36  combinations  in  which  the  nine  images  could  be  combined 
were  plotted  to  produce  a  histogram  of  matches  and  ranges,  in  which  the  effect  of 
a  few  incorrect  matches  were  typically  outweighed  by  the  true  peaks. 

The  navigation  algorithm  for  the  cart  depended  solely  on  the  vision 
system  for  it’s  inputs.  After  the  cart  took  it's  first  view  of  the  local  environment, 
it  would  move  forward  about  one  meter  and  take  another  look.  The  difference 
between  the  two  views  of  the  position  of  fixed  objects  was  the  actual  movement  of 
the  cart.  The  vehicle  would  then  plan  the  next  move,  move  one  meter,  and  then 
look  again.  Each  movement  cycle  took  approximately  fifteen  minutes. 


Problems  with  the  system  limited  it's  effectiveness.  Firstly,  the  vision 
system  could  not  reliably  identify  features  in  a  relatively  bland  environment  of 
primarily  simple  polygons.  Obstacles  which  lacked  sufficiently  high  contrast  were 
difficult  for  the  correlation  operator  to  accurately  localize.  This  had  the  effect  of 
introducing  significant  noise  into  the  system.  Secondly,  the  self  navigation  module 
was  very  sensitive  to  the  noise  effect.  The  cart  motion  was  solved  for  by  finding  a 
3D  translation  and  rotation  which  best  fit  the  before  and  after  scenes.  The  small 
errors  in  position  between  the  two  sets  of  images  would  force  the  system  to 
converge  on  an  improper  transform,  introducing  a  navigational  error  into  the 
system.  After  a  positional  error  was  introduced  into  the  system,  the  path  planner 
had  inaccurate  data  on  which  to  act  and  ultimately,  the  system  would  get  lost 
and  run  into  an  obstacle. 

2.  OSU  Hexapod 

Whereas  the  Robot  Rover  made  use  of  a  passive  triangulation  system,  the 
OSU  Hexapod  used  an  active  triangulation  system  for  vision  [Ref.  12:pp.  134-137, 
pp.  143-145].  The  purpose  of  this  vision  system  was  not  for  navigation  or  local 
terrain  identification  but,  instead,  to  designate  footholds  for  the  hexapod  to 
utilize. 

A  legged  vehicle,  as  was  discussed  before,  must  be  able  to  identify 
individual  footholds  as  it  moves  along  a  path.  Over  smooth  terrain  this  is  not 
critical.  However,  over  rough  terrain,  without  terrain  preview  information,  the 
vehicle  is  unable  to  safely  move  (except  at  very  low  speeds)  since  a  mispositioning 


of  one  of  the  footpads  could  be  catastrophic.  Lacking  an  intelligent  vision  system 
for  the  OSl*  Hexapod,  a  human  operator  had  to  designate  footholds  to  the 
hexapod.  This  was  done  by  means  of  a  laser.  The  operator  would  pick  out  likely 
footholds  within  the  range  of  motion  of  the  legs  and  illuminate  them  with  the 
laser.  Two  charge-injection-device  (CID)  cameras  were  mounted  on  masts  above 
the  vehicle  with  fields  of  view  looking  forward  and  down.  Special  hardware  was 
built  to  determine  to  position  of  the  designated  point  in  both  cameras  fields  of 
view  and  to  transmit  that  information  to  a  PDP-11/70  computer  to  calculate  the 
range  and  bearing. 

In  an  indoor  environment,  this  vision  system  proved  to  be  effective  within 
the  limits  of  the  system.  Errors  in  the  position  determination  were  1.34  cm 
forward.  1.16  cm  lateral,  and  1.2  cm  vertical.  The  half-cycle  time  of  motion  was 
25  seconds,  allowing  the  operator  sufficient  time  to  designate  the  next  foothold 
and  the  system  time  to  accept  or  reject  it.  The  motion  was  normally  smooth  with 
no  noticeable  delay.  However,  when  the  system  was  exposed  to  an  environment 
with  a  large  amount  of  background  noise  in  the  form  of  ambient  light,  the  system 
had  difficulty  distinguishing  the  designated  spot.  The  laser  used  to  designate  the 
footholds  was  a  five  milliwatt  helium  neon  laser.  Apparently,  this  was  not  enough 
power  to  provide  adequate  contrast  with  the  ambient  light. 

3.  Autonomous  Land  Vehicle 

The  present  state  of  the  art  in  range  imaging  is  represented  by  the 
Autonomous  Land  Vehicle  (ALVj  [Ref.  26:pp.  1-26].  The  ALV,  one  of  the  main 


application  areas  for  DARPA’s  Strategic  Computing  program,  is  a  program  that 
combines  recent  advances  in  sensors,  vision,  reasoning  and  advanced  computer 
architectures  with  a  mobile  8-wheeled  vehicle  in  an  attempt  to  advance  and 
demonstrate  autonomous  navigation  and  decisionmaking.  The  test  site  for  the 
ALV  encompasses  nearly  every  kind  of  terrain  an  autonomous  vehicle  could  be 
envisioned  encountering.  Initial  testing  will  only  be  on  roads  but  will  eventually 
progress  to  include  open  terrain.  At  present  the  system  is  required  to  negotiate  a 
two  kilometer  asphalt  road  at  speeds  up  to  ten  kilometers  per  hour  [Ref.  26:pp. 
8,16], 

The  vision  system  for  the  ALV  is  a  combination  of  an  RCA  color  video 
CCD  matrix  TV  camera  and  the  ERIM  laser  range  scanner  [Ref.  26:pp.  19-23]. 
The  vision  system  switches  between  the  two  imaging  systems  based  on  the 
requirements  of  the  reasoning  subsystem.  The  processing  algorithm  for  the  video 
subsystem  detects  the  difference  in  color  between  the  road  surface  and  the 
shoulder  and  vegetated  areas.  The  ranging  system  works  analogously  to  the  video 
subsystem  except,  instead  of  looking  for  a  color  contrast,  the  ranging  systems 
looks  for  a  contrast  in  the  smoothness  of  the  terrain.  This  contrast  in  the 
smoothness  of  the  terrain  results  from  the  road  being  paved  and  therefore  much 
smoother  than  the  shoulder  which  includes  small  rocks,  potholes  and  ditches  [Ref. 
5:pp.4-15].  The  advantage  of  the  ranging  system  over  the  video  system  is  that  the 
ranging  system  will  work  over  significantly  greater  variations  in  lighting  and 
weather  conditions  (day.  night,  and  even  when  the  ground  is  covered  by  a  layer  of 


snow).  While  an  algorithm  measuring  smoothness  can  do  things  like  find  obstacles 
alongside  the  road,  the  algorithm,  being  inherently  non-responsive  to  intensity  or 
color,  could  not  attempt  to  follow  a  painted  line  on  the  road  (unless  reflectance 
data  is  also  used).  Also,  while  the  video  subsystem  can  generate  an  image  in 
about  1/30 of  a  second,  the  ranging  subsystem  can  only  generate  an  image 
about  every  1/2  second.  Thus  the  two  subsystems,  video  and  ranging, 
complement  each  other  to  make  a  more  effective  real  time  vision  system. 

4.  Region  Classification 

In  1985.  Olivier  [Ref.  24:pp.  12-44],  proposed  using  the  concept  of  terrain 
or  region  classification  based  on  the  output  of  a  laser  range  finder.  His  method 
was  not  implemented  in  hardware  but  was,  rather,  simulated  on  a  PDP-11/70 
computer.  In  his  approach,  he  classified  each  terrain  cell  as  one  of  7  types:  peak, 
pit.  ridge,  ravine,  flat,  safe  slope,  and  unsafe  slope  regions. 

In  his  investigation,  Olivier  first  simulated  the  input  from  a  range  scanner 
in  terms  of  a  matrix  of  elevation  data  versus  latitude  and  longitude.  Then  a 
quadratic  equation  was  fit  to  the  local  region  about  each  surface  voxel.  The 
gradient  vector  and  Hessian  matrix  at  each  such  voxel  was  then  calculated. 
Finally,  this  data  was  used  to  generate  a  terrain  classification  for  each  voxel  in  the 
elevation  map.  From  the  results  of  this  terrain  classification,  his  algorithm  would 
print  out  a  symbolic  map  of  the  local  terrain. 

It  should  be  pointed  out  here  that  Olivier's  masters  thesis  was  specifically 
concerned  with  a  navigation  algorithm  to  do  path  planning  for  an  autonomous 


vehicle.  As  such  his  development  of  the  region  classification  algorithm  was  limited 
in  scope  and  depth. 

F.  SUMMARY 

Vision  research  for  autonomous  vehicles  has  been  a  primary  focus  of  robotic 
research  for  a  number  of  years.  The  vision  research  has  progressed  from  simple 
systems  requiring  long  periods  of  time  to  generate  a  single  images  to  almost  real 
time,  very  accurate,  complex  systems  generating  large  amounts  of  data  on  a  large 
scale  scene.  For  the  navigation  problem,  range  data  offers  one  particular  data 
format  that  can  be  used  to  simplify  the  processing  requirements  of  the  on  board 
computers.  Instead  of  trying  to  derive  individual  object  analysis  and  identification 
from  the  intrinsic  image,  the  range  data  format  can  be  utilized  to  do  region 
classification  on  local  terrain  with  each  voxel  being  one  classified  region.  In 
regards  to  the  range  imaging  systems,  the  scanning  laser  range  finder  seems  to 
offer  better  performance  than  most  other  systems  developed  to  date  and  is 
presently  the  object  of  significant  research. 

However,  any  of  the  methods  for  image  analysis  discussed  in  this  chapter  are 
as  dependent  on  mathematical  concepts  as  on  hardware.  In  the  next  chanter,  a 
detailed  analysis  of  the  mathematical  concepts  behind  the  region  classification 
method  of  image  analysis,  sis  origionally  proposed  by  Olivier  and  then  further 
extended  by  the  author,  will  be  presented  since  this  scheme  is  particularly 
applicable  to  the  laser  scanner  employed  on  the  OSU  ASV. 


III.  MODELS  FOR  TERRAIN  GENERATION  AND  CLASSIFICATION 


A.  INTRODUCTION 

j  The  first  step  in  developing  any  algorithm  to  handle  a  problem  such  as  region 

classification,  is  to  define  the  limits  of  the  problem.  In  this  case,  the  purpose  of 
the  algorithm  is  to  classify  the  visible  terrain  surrounding  an  autonomous  walking 

i 

vehicle.  Therefore,  a  terrain  data  base  has  to  be  available  to  the  algorithm.  This 
could  be  a  real  terrain  scene  that  is  imaged  by  an  actual  laser  scanner  with  the 
output  being  recorded  for  future  use  or  the  data  base  could  be  an  artificial  terrain 
scene  generated  and  made  available  to  the  computer  program.  F or  the  purposes  of 
this  thesis,  it  was  decided  to  simulate  the  terrain  using  the  available  computer 
resources. 

Once  the  terrain  has  been  generated,  an  appropriate  set  of  rules  has  to  be 
developed  to  allow  the  program  to  classify  images  in  the  field  of  view.  These  rules 
will  be  derived,  in  part,  from  the  earlier  work  of  Olivier  [Ref.  23:pp.  1-6  ;  Ref. 
24:pp.  21-44]  who  developed  a  mathematical  method  for  obtaining  the 
characteristic  first  and  second  derivatives  of  individual  voxels  in  a  range  image.  In 
this  chapter,  the  mathematical  concepts  behind  his  algorithm  will  be  developed  in 
more  detail  and  then  expanded  on. 


B.  GENERATION  OF  SYNTHETIC  TERRAIN 


There  are  any  of  a  number  of  methods  to  mathematically  generate  equally 
acceptable  synthetic  terrain  for  this  simulation.  The  choice  of  methods  depends  on 
whether  a  deterministic  or  a  random  form  of  terrain  is  desired.  Deterministic 
means  that  the  prominent  terrain  feature  of  the  scene  is  predetermined  and  then 
the  terrain  is  mathematically  generated  to  produce  a  model  of  that  feature.  On 
the  other  hand,  generation  of  random  terrain  implies  the  use  of  random  number 
generators  which  produce  entirely  random  prominent  terrain  features.  A 
deterministic  model  is  useful  in  testing  the  simulation  against  specific  terrain 
types,  such  as  an  isolated  hill,  valley,  or  cliff  terrain  feature  which  is  then  not 
effected  by  the  presence  of  other  nearby  terrain  features.  A  random  terrain  model, 
however,  is  useful  in  determining  the  overall  simulation  performance  in  a  real 
world  environment  in  which  few,  if  any,  terrain  features  are  independent  of  the 
other  surrounding  terrain  features. 

1.  Generation  of  Deterministic  Terrain 

There  are  a  number  of  possible  methods  used  to  generate  deterministic 
terrain.  The  easiest  method  is  to  simply  have  the  person  running  the  simulation 
input  the  x,  y,  and  z  values  of  the  terrain  to  be  simulated  into  a  matrix  which  can 
then  be  referenced  by  the  simulation.  This  approach  is  appropriate  for  small  areas 
of  interest,  but  for  the  simulation  of  a  large  physical  area  it  rapidly  becomes 
unmanageable. 
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To  automatically  generate  a  terrain  simulation  for  a  large  physical  area,  k 
is  necessary  to  use  some  form  of  mathematical  formula  and  then  have  the 
computer  calculate  the  terrain  matrix.  The  formula  could  be.  for  example,  a 
Fourier  series,  an  appropriately  sized  polynomial  function,  or  a  cylindrical 
function.  To  use  a  Fourier  series  to  simulate  a  specific  three  dimensional  feature, 
the  user  must  derive  the  coefficients  for  each  of  the  sinusoids  used  in  the  Fourier 
series  [Ref.  27:pp.  11-46].  For  a  simple  two-dimensional  feature,  deriving  the 
Fourier  coefficients  is  relatively  straightforward.  However,  for  a  more  complicated 
three-dimensional  feature  the  derivation  of  the  coefficients  could  prove  to  be  very 
difficult.  To  use  a  polynomial  series  would  similarly  require  the  user  to  determine 
the  coefficients  of  for  each  member  of  the  series,  an  exhaustive  and  time 
consuming  process. 

The  idea  of  a  cylindrical  function  is  a  variation  of  the  three  dimensional 
polynomial  in  which  a  two  dimensional  function  is  derived  to  simulate  the  desired 
terrain  feature  and  is  then  expanded  linearly  along  the  third  dimension  or  axis 
[Ref.  28  :pp.  34-36).  This  produces  a  stair  step  effect  and  is  adequate  to  reproduce 
terrain  features  such  as  "rolling"  terrain  and  objects  such  as  linear  "ridge  lines"  or 
"valleys".  Unfortunately,  features  such  as  hills  or  depressions  can  not  be  modeled 
by  use  of  a  cylindrical  function. 

2.  Generation  of  Random  Terrain 


One  cannot  simply  use  the  random  number  generator  to  generate  random 


elevation  coordinates  for  random  terrain  because  the  results  are  simply  a  noisy 
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but  flat  terrain  with  no  prominent  terrain  features.  Instead,  one  of  the  previously 
mentioned  techniques  is  used  with  the  exception  that  the  coefficients,  which 
determine  the  behavior  of  the  function,  are  chosen  at  random.  If  the  form  of  a 
polynomial  function,  of  sufficiently  low  order  to  keep  the  computational  time 
required  to  a  minimum,  is  chosen,  then  the  result  will  be  a  terrain  scene  ■which 
over  a  large  physical  area  will  not  be  very  realistic.  For  example,  a  function  of 
order  two  will  always  be  a  parabola  which  can  be  a  reasonable  local  simulation  of 
a  hill  or  a  depression  but  fails  over  a  large  area  because  the  resulting  parabola  will 
be  infinitely  tall  or  deep,  unlike  real  terrain  features. 

A  reasonable  simulation  is  achieved  by  use  of  a  function  similar  to  a 
Fourier  series.  In  a  true  Fourier  series,  the  sinusoids  are  based  on  the  fundamental 
frequency  and  it's  harmonics.  This  produces  a  terrain  which  is  regular  and 
repetitive.  For  example,  if  only  the  fundamental  frequency  is  used,  the  terrain 
simulation  appears  to  be  a  series  of  waves  emanating  from  one  central  point,  just 
like  the  waves  a  water  drop  makes  when  dropped  on  a  calm  surface.  Again,  the 
results  are  usable  but  not  very  realistic.  A  variation  of  the  Fourier  series  which 
produces  more  realistic  results  is  the  use  of  a  series  of  sinusoids  which,  however, 
have  random  origins,  amplitudes,  and  frequencies.  This  produces  highly  variable 
terrain  which  can  simulate  many  normal  terrain  features.  The  added  advantage  of 
using  trigonometric  functions  in  generating  the  terrain  scenes  is  that,  unlike 
orthonormal  polynomials  which  can  also  accurately  simulate  physical  terrain,  the 
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trigonometric  functions  have  properties  which  are  well  known  and  are  installed  in 
libraries  on  most  computer  systems. 

3.  Methods  Used  in  Simulation 

For  the  purposes  of  this  simulation,  it  was  decided  to  use  two  of  the 
previously  discussed  methods  to  generate  the  test  terrain  for  the  simulation.  As 
the  algorithm  was  developed,  it  was  tested  against  a  deterministic  model.  This 
was  a  local  neighborhood  test  only  (8  nearest  neighbors)  and  so  a  large  number  of 
terrain  models  could  be  manually  input  as  data.  The  algorithm  could  then  be 
validated  against  a  number  of  well  defined  terrain  types.  When  the  full  simulation 
is  run,  though,  the  area  to  be  analyzed  is  a  30  distance  unit  by  30  distance  unit 
area  (30x30  matrix)  which  is  clearly  too  large  to  be  input  manually.  Therefore  the 
simulation  uses  the  variation  of  the  Fourier  series  in  which  the  parameters  are 
randomized.  This  produces  a  large  variation  in  the  terrain  within  each  scene  and 
is  therefore  a  more  realistic  model. 

C.  QUADRATIC  APPROXIMATION  OF  SURFACE  BY  LEAST  SQUARES 
FIT 

In  his  thesis,  Olivier  bases  his  work  on  first  fitting  a  least  squares  quadratic 
surface  to  each  surface  voxel  and  it’s  eight  nearest  neighbors.  Thus,  around  each 
such  voxel,  the  elevation  (z)  of  each  point  in  the  area  of  interest  is  a  function  of 


the  x  and  y  coordinates: 
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(For  this  derivation  a  standard  right  handed  coordinate  system  is  used.  The  x 
coordinate  is  positive  to  the  right,  v  is  positive  in  a  -90°  direction  measured  in 
the  counter-clockwise  direction  and  z  (elevation)  is  positive  up.)  A  quadratic  fit  as 
a  function  of  x  and  y  has  the  following  form: 


*  2  o 

f{x,y)  =  k j  -  k2x  -  k3y  +  k^x  -  k.xy  -r  k6y‘ 


(3.2) 


The  square  of  the  error  between  the  fitted  equation  and  the  real  terrain  is  then 
defined  as 


*2  =  ES(*i  ~  k2x  +  *3 y  -  l:4x  +  ksxv  +  ke y2  ~  f(x^y)f 
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The  partial  of  the  squared  error  with  respect  to  each  constant,  k  .can  then  be 
taken  as: 
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To  minimize  the  squared  error,  each  of  the  above  equations  is  expanded,  the 
factor  of  two  is  divided  out  and  then  the  partials  are  set  to  zero.  Thus. 
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In  order  to  simplify  these  equations  so  that  they  may  be  more  easily  adapted  for 
used  by  the  algorithm  an  assumption  has  to  be  made.  This  is  based  on  the  form 
in  which  the  coordinate  data  is  formatted.  The  coordinates  of  each  voxel  in  the 
3x3  window  being  modeled  are  measured  relative  to  the  center  voxel  (the  voxel  of 
interest).  Using  this  format  the  coordinates  of  each  of  the  voxel's  is  simply  =  one 
distance  unit  from  the  voxel  of  interest.  This  has  the  added  benefit  of  non- 
dimensionalizing  the  entire  algorithm  so  that  it  can  be  applied  to  scenes  of  any 
physical  size.  TABLE  1  shows  the  resultant  relative  x  and  y  coordinates  of  the 
3x3  window'.  Using  TABLE  1.  exhaustive  calculation  shows  that  the  following 
identities  are  true  for  n$3  : 
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In  other  words,  any  terms  in  Equations  3.10  -  3.15  which  contain  a  summation  of 


an  x  or  y  coordinate  wdth  an  odd  power  will  identically  go  to  zero. 

TABLE  1 
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Accordingly.  Equations  3.10  -  3.15  can  be  rewritten  as: 
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By  again  referring  to  TABLE  1.  it  is  found  that: 

EE*2  =  6  =  EE*2 
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EE*4  =  6  =  EEv4 
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EE  *V  =  4 

*  y 
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Combining  Equations  3.22  and  3.27.  k2  is  found  to  be: 
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Combining  Equations  3.23  and  3.27,  k3  is  found  to  be: 
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And  finally,  combining  Equations  3.25  and  3.29,  kh  is  found  to  be: 
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Furthermore,  combining  Equations  3.21,  3.24,  and  3.26  with  Equations  3.27  -  3.30 
and  then  rearranging  the  resulting  equations  slightly  by  adding  the  terms 
containing  f(x,y )  to  both  sides  of  the  equations  produces  the  following  three 
equations,  respectively: 
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Dividing  Equation  3.34  by  3  and  Equations  3.35  and  3.36  by  2  produces  the 


following  equations: 
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For  the  purposes  of  this  algorithm,  finding  the  value  of  the  constant  k  is  not 
important.  (However,  to  validate  the  formula's  derived  in  this  section,  it  is 
essential  to  calculate  the  value  of  k r  This  leads  to  an  alternate  derivation,  using 
matrix  algebra  methods,  of  a  matrix  equation  to  calculate  all  the  coefficients  kl  to 
fc6.  The  method  presented  in  this  section  reproduces  the  results  obtained  by 
Olivier.  The  alternate  derivation  is  presented  in  Appendix  A.  The  section  on  the 
validation  of  the  formula’s  is  presented  later  in  this  chapter.)  Accordingly.  k.  can 
be  eliminated  from  the  above  equations  by  subtracting  Equation  3.37  from 
Equations  3.38  and  3.39  respectively.  Subtracting  Equation  3.37  from  Equation 
3.38  yields: 
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Subtracting  Equation  3.37  from  Equation  3.39  yields: 
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So  far,  this  analysis  has  yielded  5  easily  solvable  independent  equations  (reiterated 
below  as  Equations  3.31-33  and  3.40-41)  which  will  generate  5  of  the  6  constants 
necessary  to  fit  a  quadratic  equation  to  a  3x3  window  of  a  range  image  based  on  a 
least  squares  error  analysis:  the  6th  constant.  k{.  is  deemed  inconsequential  to  the 


algorithm. 
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The  purpose  of  this  algorithm  is  not  to  model  the  terrain  window  but 
rather  to  analyze  the  terrain  window  based  on  easily  definable  characteristics  of 
the  window.  These  characteristics  win  be  the  slope  and  the  curvature  of  the 
function  fitted  to  the  window.  These  characteristics  will  be  generated  using  the 
Gradient  and  the  Hessian  matrix  of  the  function  fitted  to  the  terrain  window. 
Referring  to  Equation  3.2  the  following  quantities  can  be  found.  The  first  partial 
derivative  of  f{x,y )  with  respect  to  x  is: 
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and  the  first  partial  derivative  of  f(x,y)  with  respect  to  v  is: 
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The  second  partial  of  f{x,y)  with  respect  to  x  and  v  is: 
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The  second  partial  of  /(x,y)  with  respect  to  y  is: 
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and  the  second  partial  of  /(x,y)  with  respect  to  x  is: 
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The  gradient  of  a  function  is  defined  to  be  [Ref.  29:pp.  597 
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w(*»y)  = - *  * - j 

3x  3  y 


(3.47) 


which  from  Equations  3.42  and  3.43  is: 


=  [k2  -  2  V  *  M  )  r  *  (  *3  ~  M  *  2k*V  )  3 


(3.48) 


Referring  again  to  TABLE  1.  the  gradient  of  the  voxel  of  interest  is  at  the  point 
x=0  and  y=0.  This  simplifies  Equation  3.48  to  be: 

7/(x,y)  =  k2T  -  k3j  (3.49) 
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Moreover,  it  can  be  proved  [Ref.  29:p.  597]  that  the  absolute  value  of  the 
magnitude  of  the  Gradient  of  the  voxel  is  equal  to  the  maximum  slope  at  that 
point.  Therefore  the  maximum  slope  of  the  voxel  is  defined  to  be  simply: 


8l°Pem*x  =  (  kl  +  kt  )  2 

The  Hessian  matrix  is  defined  to  be  [Ref.  30:pp.  249]: 


97(*.y) 

d2f{x,y) 

dy‘ 

dxd  y 

92/(-r,y) 

92/(r,y) 

9  ydx 

9x2 

(3.50) 


(3.51) 


By  substituting  Equations  3.44  -  3.46  into  Equation  3.51,  the  Hessian  is  found  to 
be  simply: 


H  = 


F2  k6  k 5 


*5  2*4 


(3.52) 


For  reasons  that  will  be  detailed  later  in  this  thesis,  it  is  necessary  to  determine 
the  eigenvalues  of  the  Hessian  Matrix  to  properly  classify  a  terrain  voxel.  The 
eigenvalues.  AjandA2,  of  the  matrix  are  defined  such  that  [Ref.  3l:p.  3]: 


H  -  A/ 


=  0 


(3.53) 


Using  Equation  3.52,  this  means  that: 


(3.54) 


0  = 


2k,  -  A  k. 

D  O 


k.  2k,  -  A 

o  A 


or,  solving  for  the  determinant: 


(2 *6  -  A)  (2 k4  -  A)  -  k:  =  0  (3.54) 

Expanding  Equation  3.54  algebraically  and  grouping  terms  produces: 

A2  -  (-2Jfc4-2Jfc6)A  +  {2kA2k6-kl)  =  0  (3.55) 

Using  the  quadratic  equation,  the  roots  of  Equation  3.55,  the  eigenvalues,  can  be 
found  as: 


A  = 


(2 k,  -  2k,)  x  \/{-2k4-2k6y  -  4{2k42k,-kt) 


(3.56) 


For  the  purposes  of  this  discussion,  let  the  larger  magnitude  eigenvalue  be  Aj  and 
the  smaller  one  be  Xr  To  be  able  to  apply  these  results  to  a  real  terrain  elevation 
image,  it  is  necessary  to  show  that  the  quantity  under  the  radical  is  always 
positive  and  therefore  the  eigenvalues  are  non- imaginary.  This  is  true  if: 


(-2k4-2k,)2  ^  4(2k42k,-kl) 


(3.57) 


or 


Ak24  -  8 k4k,  -  4 k2,  >  l^k4k,  -  4k\ 


(3.58) 


Subtracting  the  16 k4k,  term  from  both  sides  and  dividing  through  by  4  yields: 

(3.59) 


k“  +  k “  -  2k. k,  ^  -k, 

A  o  4  6^  d 
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This  will  always  be  true  if: 

kl  -  **  -  2 k4k6  >,  0  (3.60) 

which  must  be  true  because 

k ]  -  2 k4k6  +  kl  =  (k4  -  k6y  2  0  (3.61) 

In  summary,  due  to  the  choice  of  a  quadratic  form  for  the  fit  to  the 
window  of  voxels  surrounding  the  voxel  of  interest,  and  based  on  the  window  of 
voxel's  being  symmetrically  distributed,  the  maximum  slope  at  the  voxel  and  the 
eigenvalues  of  the  voxel's  Hessian  matrix  can  be  easily  calculated  using  Equations 
3.50  and  3.56. 

D.  TERRAIN  CELL  CLASSIFICATION  BASED  ON  QUADRATIC 
APPROXIMATIONS 

The  limited  numbers  of  ways  a  two  dimensional  quadratic  surface  can 
behave  is  the  basis  for  this  algorithm  that  classifies  terrain  cells.  The  behavior  of  a 
quadratic  surface  can  be  completely  described  by  the  gradient  (slope  as  derived 
from  Equation  3.50)  and  the  eigenvalues  of  the  Hessian  matrix  (Equation  3.56). 
To  see  how  this  is  true,  the  following  derivation  again  starts  at  Equation  3.2. 
however,  this  time  the  notation  is  changed  to  comply  with  the  notation  used  in 
Reference  30.  Acknowledgement  is  given  to  Professor  McGhee  for  allowing  the  use 
of  his  notes  in  this  section  of  this  thesis.  Equation  3.2  is  rewritten  as: 


and  in  it's  matrix  form  as: 


=  a  -  6ji  -  fc2y  -  [  x  v] 


Cll  C12 


■21 


or 


-r  t 

=  a  -  b  p  ~  p  Cp 


(3.64) 


where 


p=[*  y] 


and 


(3.65) 


C  = 


C11  C12 


C  21  C22 


(3.66) 


Without  any  loss  of  generality,  the  above  equations  can  be  simplified  by  letting 
Cp  =  c,r  The  directional  derivatives  of  Equation  3.62  can  then  be  written  as: 


d: 

—  =  -  2cux  *  2cJ2y 

ax 


(3.67) 


and 


dz 

~  ^2  ~  ~eizx  ^C22^ 
dy 


(3.68) 


In  matrix  notation  the  gradient  of  z  is  written  as: 


-?z  =  b  - 1-  2Cp 


(3.69) 
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(3.70) 

(3.71) 

(3.72) 

(3.73) 

(3.74) 


Expanding  z  in  a  three-dimensional  MacLaurin  series  about  the  origin,  z  can  be 
expressed  as: 


7 


T 

— J  rr— 


'o  *  V*  (0)p  +  —  p  Hp 
2 


(3.75) 


This  is  an  exact  expansion  of  the  quadratic  equation  even  though  the  MacLaurin 
series  is  an  infinite  series.  This  is  true  since  the  third  and  higher  order  derivatives 
of  a  quadratic  equation  are  all  identically  zero  and  the  MacLaurin  series  is  an 
infinite  series  of  increasing  order  derivatives. 

Defining  a  contour  line  as  a  line  connecting  points  of  equal  altitude  (z),  the 
equation  of  a  contour  line  can  be  derived.  For  a  given  contour  line  n,  let 
z  -  nAz  =  zn.  Then,  from  Equation  3.75.  the  equation  for  contour  line  n  is  thus: 


i  "TV-*-™ 


Note  that  without  loss  of  generality,  the  origin  of  the  x.  y  coordinate  system  can 
be  chosen  such  that  \z  =  0  (providing  that  such  a  point  exists).  If  this  is  done. 
Equation  3.76  simplifies  to 


PT Hp  =  z0  ~  -n  =  ^ 


(3.77) 


This  equation  can  be  further  simplified  by  performing  a  coordinate  rotation  in 
order  to  eliminate  any  x y  cross  product  terms.  This  results  in  a  diagonalization  of 
the  Hessian  matrix.  This  is  accomplished  by  means  of  a  similarity  transformation. 
Specifically,  suppose  that: 


X 

u 

= 

=  T 

.  y. 

.  V 

=  Tq  (3-78) 

and  T  is  defined  to  be  an  orthonormal  (rotation)  matrix.  Then,  evidently: 


T  T  T  T 

p  Hp  =  q  T1  HTq  =  q 1  \q 


where  the  middle  term  is  the  similarity  transformation  and  by  definition: 


(3.79) 


A  = 


A,  0 


0  A 


2  J 


=  TTHT  =  T  lHT 


(3.80) 


Once  again  as  described  in  the  previous  section.  Aj  and  A,  are  the  eigenvalues  of 
the  Hessian  matrix  (H)  and  T  is  called  the  matrix  of  eigenvectors  of  the  Hessian. 


The  eigenvalues  can  be  solved  for,  exactly  as  previously  described,  and  by  doing 
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so.  the  equation  for  a  contour  line  can  now  be  written  as: 


pT Hp  =  [  u  v  j 


=  A, u‘  -  X„v~  =  c 

1  ~  n 


(3.81) 


This  solution  is  based  only  on  the  surface  being  modeled  by  a  quadratic  equation, 
as  was  done  in  the  previous  section  of  this  chapter,  with  an  appropriate  choice  of 
coordinates  to  remove  the  Gradient  term  and  any  cross  products  between  x  and 
y.  The  behavior  of  Equation  (3.81)  is  dependent  exclusively  on  the  eigenvalues.  Aj 
and  A,  and  therefore  the  behavior  of  the  surface  is  a  function  of  Az  and  A0  only.  In 
examining  the  effect  of  the  eigenvalues,  it  is  seen  that  there  are  three  major 
situations  that  can  arise.  These  are  the  case  where  Aj  and  A2  are  both  of  the  same 
sign,  when  they  are  of  different  signs,  and  when  one  or  both  of  the  eigenvalues  is 
zero.  Each  of  these  cases  is  treated,  in  turn,  in  the  following  analysis. 

1.  Case  I 

Case  I  is  where  A}  and  A9  are  both  of  the  same  sign.  Here  it  is 
sufficient  to  consider  both  eigenvalues  to  be  positive  to  determine  the  shape  of  the 
contour  lines.  If  Aj  =  A2  then  Equation  3.81  becomes: 


Ai  °2  *  1,2 )  = 


(3.82) 


2  2  "  2 
u  -  f  =  -  =  r 


(3.83) 


which  is  the  equation  of  a  circle.  Furthermore,  if  the  eigenvalues  are  not  equal. 


Aj  —  A,,  then  Equation  3.81  becomes  the  equation  of  an  ellipse.  As  a  side  note, 
the  major  axis  of  the  ellipse  has  a  length  of: 


=  2 


l  A, 


(3.84) 


and  the  minor  axis  has  a  length  of: 


L  =  2 


l  A,  ! 


(3.85) 


The  sign  of  the  eigenvalues  has  the  effect  of  determining  whether  the  contour  lines 
are  growing  upward  or  downward.  If  they  are  both  negative,  then,  from  Equation 
3.81.  as  u  or  r  expand  out  from  the  origin  (get  bigger),  the  contours.  cn  becomes 
more  negative,  ie..  the  elevation  decreases,  and  so  the  point  can  be  considered  a 
hilltop.  Inversely,  if  both  eigenvalues  are  positive,  as  u  and  v  get  farther  from  the 
origin,  the  elevation  increases  and  the  point  is  a  depression. 

2.  Case  II 


Case  II  is  where  the  two  eigenvalues  have  opposite  signs.  Consider 
the  case  where  the  largest  eigenvalue.  Aj,  is  positive.  Then  the  equation  of  a 
contour  line  .  Equation  3.81.  becomes: 


Aj  a"  - 


A, 


(3.86) 


This  is  the  equation  of  a  hyperbola.  As  the  magnitude  of  the  u  coordinate 


increases  the  elevation  increases,  however,  as  the  r  coordinate  increases  the 
elevation  decreases.  Thinking  in  terms  of  a  three  dimensional  surface,  the  terrain 
analogy  to  this  surface  is  a  pass  which  is  a  low'  depression  through  a  high  terrain 
structure.  If  the  sign  of  the  eigenvalues  is  reversed,  the  increase  in  the  elevation 
due  to  the  u  term  is  less  significant  than  the  decrease  in  elevation  due  to  the  v 
term.  To  distinguish  between  the  two  situations  of  this  case,  this  situation  can  be 
thought  of  describing  a  saddle  terrain  structure.  The  difference  between  these  two 
cases  is  obtuse  but  very  real  and  needs  to  be  considered  at  length. 

3.  Case  III 

The  third  and  final  case  arises  when  one  or  both  of  the  eigenvalues 
is  zero.  Clearly,  it  is  an  impossible  situation  to  have  the  larger  eigenvalue  in 
magnitude.  Ap  equal  to  zero  and  then  to  have  the  smaller  eigenvalue  in 
magnitude.  A?,  either  negative  or  positive  in  value.  Similarly,  the  situation  when 
both  eigenvalues  are  zero  is  equally  simple  to  analyze  in  that  the  equation  of  a 
contour  line  is  identically  equal  to  zero  in  all  directions  and  therefore  the  region  is 
a  flat  or  planar  region.  This  is  not  to  be  confused  with  a  level  region.  A  flat  region 
has  little  or  no  curvature  whereas  a  level  region  has  little  or  no  slope. 

Consider  the  situation  in  which  the  smaller  eigenvalue.  A2.  is  equal 
to  zero  and  Aj  is  positive.  Referring  again  to  the  equation  of  a  contour  line,  in  the 
r  direction  there  is  no  change  in  the  elevation  while  in  the  ±u  direction  the 
elevation  increases.  This  is  a  valley  terrain  element.  Similarly,  if  Ar  the  larger 
eigenvalue,  is  negative,  then  the  terrain  remains  level  in  the  v  direction  but 


decreases  in  the  =u  direction.  This  is  the  inverse  of  a  valley  and  is  called  a  ridgt 


terrain  element.  (As  a  note,  since  nothing  in  nature  is  ever  absolutely  zero,  there 
has  to  be  a  threshold  applied  when  determining  that  one  or  both  of  the 
eigenvalues  of  the  surface  is  zero.  It  will  be  seen  later  in  the  chapter  describing  the 
results  of  the  simulation  that  this  threshold,  appropriately  called  in  the  computer 
simulation  "curv_threshold",  has  a  significant  effect  on  the  performance  of  the 
simulation.) 

TABLE  2 

CLASSIFICATION  OF  SURFACE  VOXELS 
BY  EIGENVALUES  OF  HESSIAN  MATRIX 


Aj  =  eigenvalue  with  larger  magnitude 
A,  =  eigenvalue  with  smaller  magnitude 


sign  of 
A, 

sign  of 

- 

0 

+ 

+ 

saddle 

impossible 

depression 

0 

ridge 

planar 

valley 

- 

hill 

impossible 

pass 

The  above  three  cases  represent  the  only  ways  in  which  Equation 
3.81  can  possibly  behave.  By  deciding  to  model  the  voxel  by  a  quadratic,  the 
terrain  categories  into  which  a  voxel  can  be  classified  has  been  limited  to  the 
previously  described  seven  terrain  forms.  These  are  summarized  in  TABLE  2 
where  the  sign  of  the  larger  eigenvalue  defines  the  columns  and  the  value  of  the 


smaller  eigenvalue  defines  each  row.  It  is  interesting  to  note  that  the  terrain  types 
in  column  one.  defined  by  a  negative  Ar  can  all  be  considered  as  modified  ridges. 
A  saddle  is  a  ridge  which  does  not  actually  remain  flat  on  top  but  rather  has  a 
smaller  amount  of  upward  (positive)  curvature  associated  with  it.  Similarly,  a  hill 
is  a  ridge  which  does  not  remain  level  on  top  but  curves  down  in  all  directions. 
Also,  column  three  can  be  considered  as  a  set  of  modified  valley  terrain  types.  A 
depression  (or  pit)  is  a  valley  which  does  not  remain  flat  at  it’s  bottom  but  rather 
curves  upward  in  both  directions.  And,  lastly,  a  pass  can  be  thought  of  as  a  valley 
(through  some  other  elevated  terrain  feature)  which  is  associated  with  a  small 
amount  of  downward  curvature  at  its  bottom. 

4.  The  Effect  of  Slope 

Part  of  the  above  analysis  required  setting  the  coordinate  origin  at  a 
point  such  that  the  Gradient  is  zero.  From  Equation  3.50  it  is  seen  that  the 
Gradient,  and  therefore  the  maximum  slope  of  the  voxel,  is  dependent  only  on  the 
constants  fc,  and  k 3  while  the  eigenvalues  are  based  only  on  the  constants  k 4,  k. 
and  k6.  Therefore,  the  effect  of  the  eigenvalues  is  independent  of  the  slope  and  the 
above  seven  terrain  classification  types  are  not  really  dependent  on  the  slope  being 
zero.  Logically,  if  this  is  true,  the  slope  could  be  another  way,  independent  of  the 
eigenvalues,  of  classifying  terrain.  Thus,  voxel’s  with  a  slope  equal  to  zero,  within 
some  threshold  value,  will  be  called  level  terrain.  Terrain  with  a  slope  of 
something  greater  than  the  zero  threshold  value  but  less  than  an  unsafe  slope 


threshold  will  be  designated  as  safe  slope  terrain.  Anything  else  is  called'  unsafe 
slope  terrain. 

5.  Generalized  Classification  Rules 

In  this  simulation  it  was  decided  to  use  the  slope  information  in  the 
quadratic  model  as  the  primary  classifier  for  the  voxel.  This  was  due  primarily  to 
two  reasons.  First,  in  trying  to  identify  specific  topographic  features,  the  usual 
convention  is  to  classify  terrain  based  on  the  tops  of  hills  and  ridges  or  the 
bottoms  of  valleys  and  depressions.  Small  hills  on  the  sides  of  larger  hillsides  may 
or  may  not  be  significant.  The  same  can  be  said  of  small  ridge  lines  running  along 
or  up  a  hill  side  or  down  a  valley  wall.  Secondly,  when  considering  the  concept  of 
mobility,  a  robotic  vehicle  trying  to  navigate  through  an  unpreviewed  terrain 
scene  would  be  more  interested  in  following  flat  terrain  features,  such  as  stepping 
from  hilltop  to  hilltop,  than  in  trying  to  follow  some  other  specific  terrain  features 
which  might  be  crossing  some  sort  of  steep  slope.  So  the  data  is  presented  to 
highlight  these  level  terrain  features.  This  forces  the  slope  to  be  the  most 
significant  terrain  classifier.  Admittedly,  this  is  a  somewhat  arbitrary  choice  and 
therefore  it  deserves  further  study. 

In  what  'ollows.  each  voxel  is  not  only  given  a  primary  classifier, 
but  it  is  also  given  a  secondary  classifier.  The  secondary  classifier  is  the 
information  derived  from  the  eigenvalues.  These  seven  classification  categories 
contain  the  curvature  information  of  the  terrain.  Not  only  will  an  autonomous 
vehicle  have  to  try  to  utilize  level  terrain  but  it  must  also  be  able  to  tell  the 


difference  between  specific  level  terrain  categories.  Obviously,  to  a  legged  vehicle 
there  is  a  significant  difference  between  a  flat  hilltop  and  the  fiat  bottom  of  a 
depression.  This  implies  that  each  voxel  must  be  classified  into  one  of  21  possible 
categories.  Based  on  Equations  3.50  and  3.56.  this  classification  is  possible  and 
can  be  accomplished  with  a  reasonable  degree  of  computation  ease. 

E.  VALIDATION  OF  THE  CLASSIFICATION  MODEL 

Before  attempting  to  utilize  the  algorithm  to  characterize  voxels,  as  derived 
above,  in  a  large  scale  simulation,  the  concept  had  to  be  validated.  This  was  done 
in  BASIC  on  an  ATARI  520  ST  personal  computer.  BASIC  was  chosen  over  other 
languages  due  to  it's  ease  of  use  and  rapid  debugging  capability  allowing  the  code 
to  be  written  in  a  minimum  amount  of  time.  Any  of  a  number  of  other  languages 
were  available  and  could  have  just  effectively  been  used. 

The  program  used  was  a  straightforward  implementation  using  Equations 
3.31-3.33  and  Equations  3.40-3.41  to  calculate  the  constants  necessary  for 
Equations  3.50  and  3.56.  There  was  a  second  set  of  runs  done  of  this  validation 
model  using  the  second  method  of  deriving  the  coefficients.  The  shape  of  the 
resulting  curves  was  exactly  the  same,  since  the  only  difference  was  that  with  the 
first  method  the  constant  k}  had  to  be  guessed  at  and  in  the  second  set  of  runs 
the  constant  was  calculated.  The  constant  k j  only  effects  the  vertical  displacement 
of  the  window  as  a  whole  and  not  it’s  shape. 


The  input  to  the  program  was  a  data  base  of  46  different  scenes  which  were 


generated  deterministically.  There  was  no  attempt  to  make  the  data  base 
comprehensive. 

The  methodology  worked  as  expected  with  few  exceptions.  Most  terrain 
features  were  characterized  without  error.  The  one  major  terrain  type  which  was 
unexj-  rtedly  generally  characterized  in  error  was  a  ridge  or  valley  which  was 
oriented  diagonally  across  the  terrain  scene.  When  oriented  either  port  and 
starboard,  or  fore  and  aft,  all  valleys  and  ridges  were  correctly  identified.  However 
when  they  were  diagonalized,  the  routine  identified  them  as  a  saddle  or  pass 
terrain  feature.  Upon  further  investigation,  it  w’as  discovered  that  the  quadratic 
least  squares  fit  to  such  a  terrain  scene  does  not  preserve  diagonal  valleys  or 
ridges.  By  calculating  the  elevations  of  the  voxels  by  using  the  quadratic  fit  based 
on  the  calculated  constants  ( Atj  _  g)  and  then  comparing  them  to  the  data  base 
elevations,  it  became  apparent  that  the  least  squares  fit  to  a  diagonal  valley  or 
ridge  was  indeed  a  saddle  or  pass  terrain  feature. 

The  following  are  examples  are  just  5  of  the  46  test  terrains.  In  each  table, 
the  first  entry  will  be  the  actual  terrain  type,  the  second  will  be  the  array  of 
actual  elevation  values  and  the  third  will  be  the  array  of  elevations  based  on  the 
quadratic  least  squares  fit.  After  that  will  be  the  calculated  larger  eigenvalue,  AJt 
the  smaller  eigenvalue,  A2.  the  calculated  slope,  in  degrees,  and  finally  the  terrain 
type  as  characterized  by  the  program.  For  example,  in  TABLE  3.  the  terrain  is  a 

peak  with  no  slope.  The  actual  elevation  is  the  first  3x3  matrix  and  the  calculated 
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TABLE  3 


PEAK 

NO-SLOPE 

1 

1 

1 

1 

2 

1 

1 

1 

1 

0.88 

1.22 

0.88 

1.22 

1.55 

1.22 

0.88 

1.22 

0.88 

A, 

= 

-.66 

A, 

= 

-.66 

slope 

= 

0 

PEAK 

NO-SLOPE 

TABLE  4 

RIDGE 

NO-SLOPE 

1 

2 

1 

1 

2 

1 

1 

2 

1 

1 

2 

1 

1 

2 

1 

1 

2 

1 

\ 

= 

-2.0 

a2 

= 

0.0 

slope 

= 

0 

RIDGE 

NO-SLOPE 

PASS 


PASS 


slope 


TABLE  5 


NO-SLOPE 


NO-SLOPE 


2.33 

-1.66 

0 


PLANAR 


PLANAR 


slope 


TABLE  6 


SAFE-SLOPE 


SAFE-SLOPE 
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RIDGE 


TABLE  7 
NO-SLOPE 


g 

2 

1 

0.9 

1 

1 

2 

1 

1 

0.9 

1 

2 

1 

1.84 

1.32 

0.74 

1 

1.32 

1.35 

1.32 

! 

0.74 

1.32 

1.84 

R 

K 

= 

-0.62 

i 

A2 

= 

0.48 

I 

slope 

— 

0 

I 

SADDLE 

NO-SLOPE 

terrain  is  the  second  3x3  matrix.  The  eigenvalue  with  the  larger  magnitude,  Aj,  is 
-.66  as  is  the  eigenvalue  with  the  smaller  magnitude,  Ar  The  slope  is  calculated  to 
be  0.  Finally,  on  the  bottom,  the  algorithm  classifies  the  terrain  to  be  a  peak  with 
no  slope. 


F.  SUMMARY 

The  basis  for  any  valid  computer  simulation  is  a  detailed  mathematical 
model.  The  purpose  of  this  chapter  has  been  to  logically  develop  the  model  used 
in  this  simulation  of  a  terrain  classification  scheme.  First,  the  terrain  generation 
methodology  was  discussed.  A  number  of  different  techniques  were  presented 
through  which  the  terrain  the  model  is  to  classify  could  be  generated.  Any  of  the 
methods  are  equally  as  valid  though  they  vary  in  their  ease  of  use.  The  method 
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used  in  this  simulation  is  a  variation  of  a  Fourier  series  in  which  the  sinusoids 


have  random  origins,  amplitudes,  and  frequencies.  This  method  was  chosen  since 
it  converges  rapidly  to  a  reasonable  terrain  approximation  and  due  to  the  author's 
familiarity  with  the  general  technique. 

The  mathematical  derivation  of  a  least  squares  quadratic  fit  to  a  region 
was  presented.  By  use  of  a  simplifying  assumption,  that  the  sample  points  about 
the  voxel  of  interest  are  evenly  spaced,  it  is  shown  that  the  equations  defining  the 
coefficients  of  the  quadratic  become  relatively  simple.  These  coefficients  can  then 
be  used  to  calculate  the  three  characteristic  parameters  of  a  quadratic  which  allow 
it  to  be  classified:  namely,  the  two  eigenvalues  of  the  Hessian  matrix  and  the 
slope. 

Next,  a  theory  of  quadratic  surfaces  was  presented  in  which  it  was  proven 
that  the  behavior  of  the  quadratic  surface  can  be  sufficiently  described  by  the 
slope  of  the  function  and  the  eigenvalues  of  the  Hessian  matrix  to  allow  it  to  be 
classified.  This  classification  was  divided  into  a  primary  and  secondary  category 
based  respectively  on  the  slope  and  eigenvalues.  The  primary  classification  is  into 
one  of  three  subcategories;  level,  safe-slope,  and  unsafe-slope  terrain.  The 
secondary  classification  is  into  one  of  seven  subcategories;  peak,  depression,  ridge, 
valley,  planar,  pass,  and  saddle. 

Lastly,  a  validation  of  the  algorithm  was  described.  In  that  validation  test 
one  major  discovery  was  made.  That  was  that  a  diagonal  ridge  or  valley  could  be 

incorrectly  identified  as  a  pass  or  saddle  terrain  type.  This  is  due  to  the  presence 
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of  a  crossproduct  term.  xy.  which  produces  a  sufficiently  inaccurate  least  square? 
fit  to  the  surface  that  the  program  classifies  the  terrain  erroneously.  While  the 
classification  is  erroneous,  it  should  not  be  considered  as  an  error  in  the  algorithm. 
It  is  stictiy  an  effect  of  using  a  quadratic  approximation  to  the  surface  which 
could  be  negated  by  utilizing  a  higher  order  approximation  to  the  surface. 

In  the  next  chapter  the  computer  simulation  model  will  be  discussed.  The 
model  is  written  in  LISP  and  is  run  on  an  ISI  graphics  workstation.  This  choice  of 
language  and  computer  generated  a  number  of  difficulties  in  producing  the 
simulation.  However,  these  problems  were  offset  by  the  unique  capabilities  of 
each  item. 


IV.  COMPUTER  SIMULATION  MODEL 


A.  INTRODUCTION 

Normally,  a  computer  program  is  written  to  generate  numerical  data  which 
will  then  be  either  interpreted  by  a  human  or  used  as  input  to  some  sort  of 
control  system.  Few  programs  are  written  which  are  expected  to  make  deductive 
conclusions.  This  is  because  computers  were  originally  built  exclusively  for  the 
task  of  processing  large  volumes  of  complicated  data  which  the  programmer 
would  later  analyze.  The  machine  was  responsible  for  the  number  crunching  and 
the  human  partner  was  responsible  for  symbolic  analysis.  This  was  an  efficient 
division  of  responsibilities  which  prescribed  to  each  member  of  the  team  tasks  for 
which  they  w’ere  best  suited. 

An  autonomous  vehicle  does  not  have  the  luxury  of  a  man-machine 
interaction.  The  program  which  guides  an  autonomous  vehicle  must  both 
numerically  analyze  the  data  and  also  generate  the  symbolic  conclusions  from 
w’hich  different  courses  of  action  can  be  decided  on.  In  the  case  of  this  simulation, 
the  program  is  expected  to  calculate  the  slope  and  the  two  eigenvalues  of  each 
voxel  and  then,  by  an  appropriate  set  of  rules,  decide  the  type  of  terrain  that  data 
represents.  In  this  chapter,  the  structure  of  this  program  will  be  discussed  and  the 
specific  rules  by  which  the  program  evaluates  the  data  will  be  presented.  First, 
though,  since  the  structure  of  any  program  is  as  much  a  function  of  the  desired 


results  as  it  is  of  the  choice  of  language  and  machine,  a  discussion  of  LISP  and 

• 

the  ISI  Optimum  V  Workstations  is  presented. 

B.  LISP  AND  THE  ISI  OPTIMUM  V  WORKSTATION'S 
1.  Lisp 

In  the  case  of  this  program,  the  intent  was  to  produce  a  simulation 
which  incorporated  artificial  intelligence  techniques  and  also  produced  a  graphical 
rather  than  a  tabular  output.  In  the  field  of  artificial  intelligence  there  are  two 
major  languages  in  use  today.  These  are  PROLOG  and  LISP.  Each  language  has 
it's  own  particular  strengths  and  weaknesses.  The  choice  of  LISP  as  the  language 
in  which  this  simulation  would  be  written  was  based  primarily  on  two 
considerations.  First,  LISP  (LISt  Processing)  [Ref.  32:p.  xi]  is  very  strong  in  it's 
ability  to  manipulate  lists  as  well  as  perform  numerical  analysis.  Since  an  image 
can  be  thought  of  as  list  of  pixel  values,  it  appeared  that  the  list  handling 
functions  in  LISP  would  be  of  value.  Secondly,  the  ultimate  purpose  of  this 
program  is.  potentially,  to  be  incorporated  as  a  tool  in  a  larger  program  to  do 
path  planning  and  navigation  for  the  ASV  (adaptive  suspension  vehicle)  [Ref. 
33:pp.  66-71].  A  proposed  computer  upgrade  of  the  the  ASV's  systems  includes  a 
high  performance  LISP  computer  for  this  navigation  module.  Therefore.  LISP 
appeared  to  be  the  logical  choice  of  languages  in  which  to  write  this  simulation. 
However,  Franz  LISP  (the  dialect  of  LISP  installed  on  the  ISI  Optimum  V 
Workstations  at  the  Naval  Postgraduate  School)  has  no  inherent  graphics 


capability.  This  was  a  major  concern  with  respect  to  the  desire  for  a  graphical 
output  from  the  simulation. 

As  it  turned  out.  the  choice  of  LISP  was  a  very  fortunate  one  due  to 
the  extreme  versatility  of  the  language.  To  justify  that  statement  it  is  necessary  to 
first  explain  a  few  of  the  details  of  how  LISP  operates  [Ref.  32:pp.  1-385).  LISP  is 
usually  executed  as  an  interperted  language  and  is  probably  best  described,  as  a 
language,  as  a  set  of  list  manipulating  functions.  LISP  functions  are  written  in 
terms  of  s- expressions  which  stands  for  symbolic  expressions  [Ref.  32:p.  2).  An  s- 
expression  is  generally  some  sort  of  command  or  statement  which  is  surrounded 
by  a  set  of  parentheses.  (Symbolic  expressions  are  a  form  of  a  list  wrhere  a  list  is 
defined  as  a  sequence  of  objects  inside  a  pair  of  parentheses.)  For  example, 
"(  times  8  3  )"  is  a  LISP  s-expression  which  tells  the  lisp  interperter  to  multiply  8 
by  3.  As  the  LISP  interperter  reads  an  s-expression  it  goes  into  an  evaluation 
process  where  it  attempts  to  carry  out  the  process  requested  in  the  s-expression 
through  an  evaluation  of  each  item  in  the  list  which  makes  up  that  s-expression. 
The  arguments  of  the  s-expression  (the  elements  of  the  list  which  follow  the  first 
item)  can  be  either  an  atom  or  another  s-expression.  An  atom  is  something  which 
evaluates  to  itself  (eg.,  the  number  8  in  contrast  to  another  list  or  s-expression).  If 
the  arguments  are  s-expressions  then  LISP  identifies  them  as  such  and  attempts 
to  evaluate  them.  Thus  LISP  proceeds  down  through  an  s-expression  attempting 
to  evaluate  each  embedded  s-expression  until  it  gets  to  the  bottom  level  where  all 
the  arguments  to  the  s-expressions  are  atoms.  It  then  returns  the  value  of  the 


bottom  s-expression  up  to  the  next  lower  level  of  s-expression  so  that  that  level 
can  be  evaluated.  Progressing  from  the  bottom  level  up  to  the  top  level  of  the  s- 
expression.  LISP  evaluates  each  level  of  s-expression  based  on  the  value  returned 
by  evaluating  the  next  lower  level  until  it  finally  evaluates  the  top  level  and 
returns  that  value  [Ref.  32:pp.  1-14]. 

LISP  functions  are  invoked  in  terms  of  s-expressions.  Referring 
again  to  the  previous  example  of  an  s-expression  which  is  a  mathematical 
function,  namely,  "(  times  8  3  )",  LISP  attempts  to  evaluate  the  s-expression  by 
evaluating  the  the  list.  It  assumes  that  .the  first  item  in  the  list  is  a  function  name. 
To  perform  the  requested  function.  LISP  then  attempts  to  evaluate  the  second 
and  subsequent  s-expressions.  To  evaluate  those  items  it  tries  to  match  them  with 
similar  symbols  on  what  is  called  the  oblist  which  stands  for  object  list.  If  one  of 
those  symbols  is  not  part  of  the  oblist  then  LISP  doesn’t  known  what  to  do  with 
it  and  it  returns  an  error  such  as 

Error:  Unbound  Variable:  (symbol  name) 

saying  that  the  item  is  not  bound  to  any  particular  item  on  the  oblist.  When  a 
symbol  is  bound,  it  is  interned  on  the  oblist  [Ref.  32:p.  236].  This  is  done  by 
adding  the  symbol  to  the  oblist  and  associating  it  with  a  pointer  which  points  to 
the  address  in  memory  of  either  a  routine,  list,  or  constant  representing  the  value 
of  that  particular  symbol.  Returning  to  the  example,  as  the  function,  in  the  form 


of  an  s-expression.  is  evaluated  by  LISP.  LISP  finds  the  pointer  on  the  oblist 
which  corresponds  to  "times".  LISP  then  jumps  to  the  routine  in  memory  at  that 
pointer.  This  routine  takes  the  next  two  items  in  the  s-expression.  adds  them  and 
then  returns  the  sum  as  the  value  of  that  s-expression. 

What  makes  LISP  such  a  dynamic  language  is  that  the  oblist  is  not 
a  static  structure.  All  variables  as  they  are  input,  all  lists  which  are  created,  all 
strings,  etc.  have  to  be  accessible  through  the  oblist  and  so  the  oblist  is  a 
structure  which  is  constantly  in  flux.  Not  only  can  lists  and  atoms,  for  example, 
be  added  to  the  oblist,  but  so  can  a  user  defined  function,  since  a  function  is 
actually  nothing  more  than  a  list  of  calls  to  other  LISP  functions.  In  fact,  most 
LISP  programming  is  done  by  a  user  writting  his  own  list  of  calls  to  the  functions 
which  perform  the  sequence  of  manipulations  he  requires.  Then  that  list  is 
interned  on  the  oblist  under  the  name  the  user  wishes  to  call  his  new  function. 
This  is  done  through  the  LISP  function  defun  which  stands  for  define  function. 
After  that  the  user  can  then  invoke  this  new  function  in  exactly  the  same  manner 
as  any  other  LISP  function  (by  writing  it  in  an  s-expression  along  with  it’s 
arguments  and  then  evaluating  the  s-expression).  (This  process  of  defining 
functions  is  more  completely  described  in  Reference  32.  Chapter  3.) 

Not  only  can  the  programmer  write  a  new  function  in  terms  of  calls 
to  other  LISP  functions,  but  a  new  function  can  also  be  created  by  calls  to  foreign 
functions  .  i.e.,  functions  written  and  compiled  in  "Fortran".  "C",  or  "Pascal". 
This  capability  contributes  substantially  to  LISP's  great  versatility.  It  is 


implemented  through  a  LISP  function  called  cfasl.  To  utilize  this  capability,  the 
foreign  function  has  to  be  written  in  one  of  the  previously  mention  languages  as 
subroutine  or  function  (eg.,  without  an  "main"  routine)  and  then  compiled  into 
object  code,  without  being  assembled,  by  invoking  the  compiler  with  appropriate 
options  (  -c  for  the  "Fortran"  and  "C"  compilers  in  Unix)  [Ref.  34:pp.  8-4  to  8-8]. 
When  the  "cfasl"  function  is  called,  LISP  reads  the  object  code,  assembles  it,  and 
then  puts  the  pointer  to  the  function  into  the  oblist,  bound  to  the  name  the  user 
wishes  to  use  to  call  the  new  function.  From  that  point  on  the  "foreign"  function 
appears  be  a  LISP  function.  The  fact  that  it  is  written  in  another  language  and 
can  access  the  unique  libraries  of  that  language  is  totally  transparent. 

This  capability  is  not  well  documented  in  most  LISP  reference  books 
and  there  are  a  few  peculiarities  involved  in  it’s  use.  First,  arguments  to  "foreign" 
functions  are  passed  as  pointers.  This  usually  causes  no  problem  and  the 
programmer  is  often  unaware  of  the  way  these  arguments  are  passed.  In  Fortran 
variables  are  passed  as  pointers  anyway,  so  no  special  efforts  must  be  taken  when 
"cfasl"ing  Fortran  routines.  In  "C",  though,  there  is  no  set  convention  for  passing 
arguments.  Arguments  can  be  passed  by  passing  either  the  item  or  by  passing  a 
pointer  to  the  item.  A  problem  specifically  arises,  in  relation  to  the  graphics 
routines,  when  attempting  to  pass  the  Bitmap  Descriptor.  BMD  ,  which  is  a 
pointer  to  the  space  in  memory  containing  a  bitmap.  LISP  passes  this  parameter, 
when  cfasl’ing,  by  passing  a  pointer  to  a  memory  location  where  the  BMD.  a 
pointer  itself,  is  stored.  A  "C"  program  attempting  to  utilize  the  BMD  passed 


from  the  LISP  environment  has  to  declare  the  argument  it  receives  as  a  pointer. 
When  the  "C"  program  attempts  to  retrieve  the  BMD  by  accessing  the  data  in 
the  memory  location  it  must  also  declare  the  new  item  as  a  pointer  type  variable. 
This  forces  the  "C"  compiler  to  generate  a  warning  since,  apparently.  "C" 
compilers  do  not  like  pointers  to  pointers.  The  warning  will  not  abort  the 
compilation.  However,  the  compiler  will  protest  vehemently. 

Secondly.  LISP  stores  "characters"  and  "strings"  of  characters 
differently  than  it  does  an  atom.  Consequently,  when  passing  a  "string"  or 
"character"  as  an  argument  to  a  foreign  routine  (or  a  LISP  function  for  that 
matter),  the  programmer  must  be  careful  to  not  attempt  to  pass  the  atom 
representing  those  variables  instead  of  the  "string"  itself.  For  example,  if  the 
variable  "title"  in  LISP  has  been  bound  to  the  string  ’NOISY,  then  there  is  a 
difference  between  the  following  two  LISP  calls  to  a  foreign  graphics  function 
which  writes  a  string  to  a  bitmap: 

(paint  character  string  "NOISY") 

and 

(paint  character  string  title) 

The  first  s-expression  will  work  since,  in  LISP,  anything  inside  double  quotes  is  a 
"string".  However,  the  second  will  generate  an  error  and  will  not  work  since  the 


variable,  title,  is  stored  as  -an  atom.  To  pass  a  variable  "string"  to  a  foreign 
function,  the  programmer  can  make  use  of  the  LISP  function  which  converts 
atoms  into  "strings".  In  FRANZ  LISP  this  function  is  "(getpname 
variable_name)"  [Ref.  32:p.  251]. 

In  summary.  LISP  is  a  very  versatile  and  rich  language.  For  image 
processing,  its  list  processing  functions  provide  a  good  set  of  resources  to  utilize. 
For  mathematical  manipulations,  it  has  a  good  set  of  tools  with  which  to  work.  It 
can  access  "foreign"  functions  more  easily  than  most  other  languages  and  passing 
parameters  to  these  functions  is  direct  and,  for  the  most  part,  easily  accomplished. 
Finally,  though  it  is  a  slow  language  in  it’s  interperted  form  on  most  "normal" 
machines,  there  are  a  number  of  new  "LISP"  machines  which  promise 
performance  equal  and  possibly  superior  to  the  performance  of  other  compiled 
languages  on  present  day  high  speed  computers. 

2.  ISI  Optimum  V  Workstations 

The  computer  on  which  it  was  decided  to  implement  this  program  is 
an  ISI  Optimum  V  Workstation.  This  is  a  graphics  workstation  with  a  Motorola 
68010  microprocessor  (since  upgraded  to  a  68020  microprocessor),  two  megabytes 
of  available  memory  and  which  htis  the  UNIX  4.2BSD  operating,  system  installed. 

The  workstations  feature  a  high-resolution  bitmapped  display  for 
graphic  display.  The  addressable  display  area  on  the  screen  is  1280  pixels  width 
by  1024  pixels  height  on  a  19  inch  tube.  As  implemented  on  the  available 
workstations,  the  display  is  Black  and  White,  although  color  workstations  are  also 


available.  In  order  to  use  the  graphics  environment  available  on  the  workstations, 
a  window  system  allowing  the  presentation  of  multiple  overlapping  windows  has 
been  incorporated  into  the  4.2  BSD  UNIX  operating  system.  Each  window  can  be 
configured  as  one  of  a  number  of  a  standard  terminal  types  with  the  default 
window  type  being  an  ANSI  x3.64  terminal  running  the  standard  4.2  BSD  UNIX 
tty  line  discipline.  The  graphics  interfaces  included  in  these  machines  are' 
SIGGRAPH  Core  and  GKS  (Graphical  Kernel  System)  which  are  two  of  the 
standardized  graphics  interfaces  available  with  considerable  device  independence 
[Ref  ,34:pp.  1-1  to  1-8]. 

Besides  the  fact  that  these  systems  were  readily  available  for  this 
research,  the  reason  for  deciding  on  their  use  is  their  outstanding  graphics 
capabilities.  Installed  as  libraries,  which  are  easily  accessible  to  the  user  through 
either  Fortran  and  C,  are  a  standard  set  of  tools  allowing  the  user  to  perform  a 
wide  variety  of  graphical  manipulations  [Ref.  35:p.  iii].  The  libraries  are  divided 
into  three  levels  depending  on  the  user’s  requirements.  The  upper  level  is  the 
libtooh  library.  This  library  provides  to  the  programmer  a  set  of  easy  to  use 
routines  as  interfaces  for  graphics  programs.  Typically,  included  at  this  level  are 
routines  to  open  and  close  new  windows,  create  and  manipulate  menus,  and  to 
change  cursor  patterns  [Ref.  36:pp.  2-1  to  2-6).  With  the  exception  of  opening  and 
closing  windows,  little  use  was  made  of  this  library.  The  next  level  of  library  is 
the  libvt.  which  stands  for  library  virtual  terminal.  Each  open  window  is  a ' 


different  terminal  with  it's  own  set  of  operating  characteristics,  even  though  it  is 


just  one  part  of  the  display  on  the  physical  terminal  (the  desktop  ).  Accordingly, 
each  window  is  termed  a  virtual  terminal.  The  libvt  provides  the  programmer 
with  a  number  of  primitives  for  window  management,  pointer  positioning, 
attribute  management,  object  painting,  clipping  management,  region 
manipulation,  refresh,  cursor  management,  font  management,  input  buffering,  and 
output  buffering  [Ref.  37:p.  l-l].  With  these  primitives,  the  programmer  is 
provided  an  extensive  library  to  generate  the  required  graphical  output  on  a 
virtual  terminal.  At  this  level,  though,  the  programmer  does  not  have  access  to 
the  actual  bitmap  which  is  being  displayed  on  the  virtual  terminal.  He  can  only 
invoke  the  specific  library  routines  provided  which  will  then  perform  the  required 
manipulations  on  the  bitmap  which  is  associated  with  a  particular  window.  This 
relieves  the  programmer  of  the  responsibility  of  having  to  manipulate  individual 
bits  on  a  display.  For  this  research,  there  existed  a  requirement  for  printed  copies 
of  the  displays  to  be  generated  for  inclusion  in  this  thesis.  At  the  libvt  level,  there 
are  no  routines  available  which  will  dump  copies  of  the  displays  into  files  which 
can  then  be  printed  on  a  printer.  This  requires  access  to  the  actual  bitmap. 

The  lowest  level  of  library  in  the  ISI  graphics  system  is  called  libbm 
which  stands  for  library  bitmap.  This  is  the  library  which  is  most  extensively  used 
in  this  simulation  to  generate  and  manage  the  graphics  required  for  the  output. 
Libbm  is  the  library  in  which  the  programmer  has  access  to  the  individual 
bitmaps.  Most  of  the  capabilities  in  the  libvt  are  present  in  libbm  with  extra 
routines  which  are  appropriate  when  working  with  bitmaps  [Ref.  38:pp.  2-1  to 
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2-16  .  A  bitmap  is  a  portion  of  memory  which  holds  the  data  the  video  hardware 
uses  to  produce  a  picture  on  the  screen.  To  produce  a  display  on  a  screen,  the 
video  hardware  turns  individual  pixels  on  by'  illuminating  the  pixel  with  an 
electron  beam  from  the  electron  gun  at  the  back  of  the  picture  tube  or  off  by  not 
illuminating  the  pixel.  The  beam  is  scanned  across  the  picture  tube  in  a  raster 
scan  starting  from  the  upper  left  corner  down  to  the  lower  right  corner,  row  by 
row.  The  video  hardware  turns  the  electron  beam  on  and  off  according  to  the 
sequence  of  bits  in  memory,  an  "on”  bit  (1)  turning  the  beam  on  and  an  "off’  bit 
(0)  turning  off  the  beam.  This  sequence  of  Ts  and  0's  in  memory  is  called  a 
bitmap  since  it  is  the  "map"  by  which  the  electron  gun,  which  draws  the  picture 
on  the  picture  tube,  is  controlled.  The  mechanics  by  which  the  libbm  library 
actually  sets  and  resets  bits  in  the  map  will  not  be  discussed  here  but  the  reader  is 
directed  to  Reference  34  for  further  details. 

The  libbm  allows  a  programmer  to  allocate  a  block  of  memory  in 

which  the  bitmap  will  be  maintained  and  to  display  that  bitmap  in  any  available 

open  window.  The  library  also  has  a  set  of  routines,  similar  to  the  libvt,  which 

allows  the  programmer  to  perform  graphics  manipulations  such  as  drawing  a  line 

from  point  A  to  point  B,  setting  the  position  of  the  cursor  pointer,  setting  the  line 

style,  settir  the  background  color,  and  more.  By  use  of  these  routines,  the 

programmer  is  relieved  of  the  responsibility  of  calculating  which  bits  need  to  be 

changed  to  perform  these  actions  and  how  these  changes  effect  each  byte  in  the 

bitmap,  yet.at  the  libbm  level,  the  programmer  does  retain  the  option  of 
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manipulating  individual  bits.  The  one  routine  which  dictated  the  use  of  the  libbm 
library  instead  of  the  libvt  library  is  called  "bmReadBit"  through  which  the 
individual  bits  in  the  bitmap  can  be  read.  This  allows  a  function  written  in  "C" 
to  sequentially  read  the  entire  bitmap  of  a  display  created  by  the  LISP 
simulation,  convert  it  into  hexadecimal  code,  and  then  write  that  code  to  a  file 
which  a  laser  printer  can  then  use  to  print  a  copy  of  the  display.  That  function, 
called  "tofile.c",  is  then  "cfasl"ed  into  a  LISP  function  which  the  simulation  uses 
to  produce  the  graphical  output.  Of  the  foreign  graphical  functions  written  for 
this  simulation,  "tofile.c"  is  the  only  one  written  in  "C". 

A  number  of  other  foreign  graphical  routines  were  written  and  then 
"cfasl"ed  into  LISP  (again  attesting  to  the  versatility  of  LISP)  to  perform  the 
graphics  required  in  this  simulation.  These  graphical  routines  were  all  written  as 
Fortran  subroutines  and  integer  functions.  These  are  included  as  Appendix  B. 
Each  routine  is  documented  in  this  appendix  as  to  how  to  use  it,  how  it  is 
invoked,  what  parameters  must  be  passed,  and  how  it  must  be  compiled. 
Consequently,  these  matters  will  not  be  further  discussed  here. 

C.  STRUCTURE  OF  THE  PROGRAM 

This  LISP  program  has  been  written  as  a  file  which  contains  the  constants, 
user  defined  functions,  etc.,  which  perform  the  required  actions.  It  is  titled  as 
LSCAN  (Laser  SCANner)  and  it  is  loaded  into  LISP  with  the  command  "(load 
’lscan)".  A  copy  of  the  file  is  included  as  Appendix  C.  Every  attempt  has  been 


made  to  adequately  document  the  purpose  of  each  function,  constant,  and  section 
and  to  keep  the  code  reasonably  readable  without  burying  the  code  in  comment. 
The  file  is  written  in  four  sections.  The  first  section  "cfasl"s  the  graphics  routines 
into  the  oblist.  The  actual  format  of  the  "cfasl"  command  is  described  in  the 
introduction  to  this  section.  The  second  section  loads  the  constants  which  are 
referenced  by  the  different  user  defined  functions.  The  purpose  of  each  constant  is 
described  as  it  is  entered.  The  third  section,  which  is  the  largest  section,  is  where 
all  the  user  defined  functions  are  written.  As  LISP  loads  this  portion  of  the  file, 
each  function  is  sequentially  interned  on  the  oblist.  Included  as  a  header  to  each 
function  is  a  short  definition  of  that  function's  purpose.  Finally,  the  section  which 
would  normally  be  called  the  main  program  is  presented.  It  consists  of  four  calls 
to  the  function  "(def  window)"  and  then  a  call  to  the  function  "(keep  running)". 
Each  time  "(def  window)"  is  called  it  defines  the  parameters  of  a  window 
(window  width  and  window  height  in  pixels,  window  position  in  terms  of  x  and  y 
coordinates  of  the  upper  left  corner  of  the  window  in  pixels  relative  to  the  position 
of  first  window,  and  the  window  title)  and  attaches  this  list  of  parameters  as 
properties  to  the  window  name."  wl"  through  "w4".  Finally  "(def  window)" 
appends  the  window  name,  along  with  it's  properties,  to  the  end  of  the  list 
"windows".  The  function  "(keep  running)"  calls  the  function  "(drawl)"  and  then. 


iteratively  calls,  the  function  "(redraw)"  until  the  appropriate  number  of  runs  of 
the  simulation  have  been  completed. 


The  functions  "(drawl)"  and  "l redraw)"  are  very  similar  except  that 
"(drawl)"  initializes  the  bitmaps  on  which  the  graphical  output  is  drawn  while 
"(redraw)"  only  has  to  clear  each  bitmap  of  the  previous  image  so  that  the  new 
one  can  be  drawn.  Figure  4  is  a  flow  chart  of  these  two  function  without  the  the 
initialization  or  bitmap  clearing  routines  included. 

The  routine  "(initialize)"  needs  to  be  explained  because  it  establishes  the 
way  the  windows  display  the  graphical  output.  Initialize  takes  the  information  in 
the  list  "windows"  and  opens  a  window  for  each  window  on  the  list.  The  way  the 
program  is  presently  written,  there  are  four  windows  open.  The  first  window  is 
"wl"  and  it  is  the  reference  for  positioning  the  other  three  windows.  It  is  titled 
"TERRAIN"  and  its  upper  left  corner  is  at  the  coordinates  stored  in  the  constants 
"wlx"  and  "wly".  The  second  window,  "w2",  is  titled  "NOISY  TERRAIN"  and 
it  is  directly  below  "wl".  The  upper  left  corner  of  "w2"  is  directly  in  line  with 
"wlx"  but  displaced  down  by  the  sum  of  the  window  height  and  the  desired 
vertical  separation  between  windows  (the  formula  in  the  second  call  to 
"(def  window)"  ).  These  constants  are  stored  as  "wwheight"  and  "ysep" 
respectively.  The  position  of  windows  "w3"  and  "w4"  are  to  the  right  of  "wl"  and 
"w2"  respectively.  They  are  titled  "CLASSIFICATION  MAP"  and  "NOISY 
CLASSIFICATION  MAP".  By  not  writing  "(initialize)"  with  four  set  of 
constants,  but  rather  having  "(initialize)"  get  the  position  of  the  four  windows 
from  the  list  "windows",  the  symbolic  output  of  the  simulation  can  be  easily 
modified.  New  window  width’s,  window  height's,  and  vertical  and  lateral  window 


Figure  4. 

Drawl  and  Redraw  Flowchart 


separation  can  be  changed  by  simply  changing  the  constants  wwidth.  wheight. 
wlx.  wly.  xsep.  ysep.  or  by  changing  the  formula's  by  which  the  four  calls  to 
(def_window).  in  the  "main  program"  section,  calculate  the  x  and  v  coordinates 
of  each  window.  This  allows  the  user  the  maximum  flexibility  in  controlling  the 
displays  through  easily  modified  constants  and  routines. 

The  function  "(draw  terrain)"  is  just  one  of  the  many  functions  which  are 
exclusively  concerned  with  the  graphical  display  of  the  data.  They  are 
straightforward  consisting  primarily  of  calls  to  set  the  cursor  position  and  to  draw 
a  line  to  another  coordinate.  These  functions  are  well  documented  and  therefore 
will  not  be  discussed  further  except  as  follows.  Included  in  the  graphics  functions 
are  a  set  of  functions  whose  purpose  it  is  to  draw  (paint)  characters  on  a  bitmap 
through  a  series  of  calls  to  a  function  to  draw  the  required  set  of  straight  lines. 
The  simulation  uses  these  characters  to  symbolically  identify  each  terrain  voxel  on 
the  map.  There  are  ten  of  these  functions  which  effectively  draw  an  ordinary  ascii 
character  to  the  bitmap.  These  functions  were  written  before  the  anomaly  about 
passing  strings  and  characters  to  foreign  functions  was  discovered.  There  w'as 
another  lisp  function  written  before  those  ten  functions  wrere  written  whose 
purpose  it  was  to  paint  a  character  to  the  bitmap  by  accessing  a  function  in 
libbm.  "bmPaintChar".  That  lisp  function  never  worked  because  an  atom  instead 
of  a  character  was  being  passed.  These  functions  could  now  easily  be  replaced  by 


The  function  "(cop>  bm  tofiiel"  utilizer  the  foreign  function,  "to  file", 
discussed  in  the  previous  section  to  convert  and  copy  the  bitmap  to  a  file  that  the 
laser  printer  can  print.  The  only  new  thing  "(copy  bm  tofile)"  does  is  to 
increment  a  counter  (named  filenum)  and  concatenate  the  counter  to  the  string 

"figure.  ".  This  creates  a  neve  file  name  into  which  the  function  will  put  the 

bitmap.  By  concatenating  filenum  and  "figure.  "  together,  they  are  converted 
into  an  atom  and  must  be  reconverted  into  a  string. 

The  function  "(make_real_terrain)"  is  where  the  artificial  terrain  is 
generated.  The  function  "(make  real  terrain)"  uses  the  random  number  generator 
to  produce  five  random  numbers.  The  random  numbers  are  rl  and  r2  which  are 
the  x  and  y  coordinates  of  the  origin  of  the  sinusoids.  r3  which  determines  the 
frequency  of  the  sinusoid,  r4  which  is  the  amplitude  of  the  sinusoids,  and  r5  which 
is  the  number  of  alternating  sines  and  cosines  that  will  be  summed  together  to 
produce  the  elevation  map.  The  summation  of  these  sinusoids  is  stored  in  the 
array  "real  terrain".  The  function  "(make_noisy_terrain)"  simply  adds  a  random 
number  between  ±  noise  limit.  a  constant  defined  in  the  first  section  of  the 
program,  to  the  elevations  in  the  array  "real  terrain"  and  then  stores  these  new 
elevations  into  the  array  "noisy  terrain"  to  impress  a  random  noise  signal  on  top 
of  the  terrain. 

The  terrain  classification  is  done  within  the  functions 
"(draw _rclass_map)".  meaning  draw  real  (no  noise)  symbolic  terrain  map.  and 
"(draw_nclass_map)",  meaning  draw  noisy  symbolic  terrain  map.  by  calls  to  the 
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functions  "( make _rc lass  map)"  and  "(makenclassmap)".  Figure  5  presents  a 
flow  chart  of  "(draw  rclass  map)".  Since  each  of  the  functions  operate  in  exactly 
the  same  manner,  only  the  function  "(draw_rclass_map)"  is  flowcharted  and 
discussed  in  this  chapter.  To  analyze  the  function  "(draw  nclass  map)".  Figure  5 
is  applicable  with  the  exception  that  the  name  "rclass"  is  replaced  by  the  name 
"nclass".  The  function  "(drawrclassmap)"  first  calls  the  function 
"(make  rclass  map)"  illustrated  by  Figure  6.  This  function  returns  a  list  called 
' relass  map  in  which  data  on  each  voxel  is  stored  as  a  set  of  properties.  These 
properties  are  attached  to  one  of  900  atoms  which  make  up  the  list  ‘rclass  map. 
These  atoms  are  named  ’pixel-1  through  ’pixel-900.  (There  are  900  members  to 
this  list  because  the  terrain  scene  that  is  being  classified  is  30  X  30  voxels  in  size.) 
One  of  the  properties  attached  to  each  pixelname  is  a  property  called  ’letter.  This 
defines  the  symbol  that  will  be  drawn  on  the  symbolic  map.  The  function 
"(draw  rclass  map)"  sequentially  checks  each  member  of  the  list  ‘rclass  map  for 
it’s  ’letter  and  then  draws  that  symbol  onto  the  bitmap  at  the  position  associated 
with  that  voxel.  Once  the  entire  list  has  been  checked  and  all  the  symbols  drawn 
on  the  bitmap,  "(draw_rclass_map)"  displays  the  bitmap  in  the  designated 
window  on  the  screen. 

The  work  of  classifying  the  terrain  voxels  is  done  in  "(make  rclass  map)" 
and  "(make  nclass  map)".  In  these  functions,  the  arrays  in  which  the  terrain 
maps  are  stored  are  read  and  the  x,  y.  and  z  coordinates  of  each  voxel  are 
attached,  as  properties,  to  the  voxel  name,  pixel-1  through  pixel-900.  As  each 


array  is  read.  "(make_rclass  terrain)"  checks  to  see  if  the  voxel  is  an  edge  voxel, 
which  is  unclassifiable  (since  it  does  not  have  a  complete  set  of  nearest  neighbors 
available).  If  it  is.  then  it  assigns  it  a  letter  of  "E",  for  edge,  and 
"(make_rclass_map)"  continues  on  to  the  next  sequential  voxel  in  the  array.  If 
the  voxel  is  not  an  edge  voxel,  "(make  rclass  map)"  invokes  the  function 
"(classify)"  and  passes  to  it  the  array  name,  "real_terrain",  and  the  voxel’s 
coordinates  in  the  array,  pi  pj.  The  function  "(classify)"  then  calculates  the  slope 
and  the  two  eigenvalues  for  that  voxel,  from  the  previously  derived  formulas,  and 
by  using  the  rules  developed  in  the  last  chapter,  classifies  the  voxel.  When  it 
classifies  the  voxel,  "(classify)"  attaches  to  the  voxel  name  the  slope,  the  two 
eigenvalues,  a  primary  classifier,  a  secondary  classifier,  and  a  the  symbolic  map.  It 
is  this  ’letter  which  the  function  (drawrclassmap)  accesses  to  decide  what  to 
paint  to  the  map. 

After  "( make _rc lass  map)"  and  " (make  nclass  map) "  have  completely 
classified  and  drawn  their  respective  symbolic  maps  to  "w3"  and  "w4", 
"(keep  running)"  checks  to  see  if  the  required  number  of  runs  of  the  simulation 
have  been  executed.  If  they  have,  the  function  calls  "(finish)"  which  closes  the 
four  windows,  deallocates  the  space  in  memory  reserved  for  the  four  bitmaps,  and 
then  exits  to  the  LISP  environment.  If  there  are  more  runs,  "(keep  running)" 
invokes  "(redraw)"  which  clears  the  old  bitmaps  and  then  starts  over  again 
making  a  new  set  of  real  and  noisy  terrain. 


D.  CLASSIFICATION  AND  SYMBOLS 

The  methodology  for  classifying  individual  pixels  was  discussed  in  the 
previous  chapter.  To  reiterate,  the  primary  classifier  of  a  voxel  is  it's  slope.  The 
slope  can  either  be  level,  safe  slope,  or  unsafe  slope.  The  thresholding  between 
each  regime  is  at  the  users  discretion  but,  as  will  be  shown  in  the  next  chapter, 
the  thresholds  have  a  major  effect  as  to  the  usability  of  the  output.  The  secondary 
classifier  is  based  on  the  value  of  the  two  eigenvalues  of  the  Hessian  matrix,  the 
larger  magnitude  eigenvalue  being  Aj  and  the  smaller  being  Ar  In  the  case  of  the 
secondary  classifier,  there  are  seven  possible  classes  and  one  threshold.  The 
threshold  is  the  limit  between  a  zero  eigenvalue  and  a  nonzero  eigenvalue  and,  in 
the  simulation,  it  is  called  "curv  threshold".  This  leads  to  a  possibility  of  21 
different  classes  into  which  each  voxel  can  be  put.  If  a  different  symbol  were  to  be 
assigned  to  each  type  class  the  display  might  be  difficult  for  a  human  to  interpret. 
Consequently,  a  scheme  had  to  be  developed  which  would  reduce  the  number  of 
symbols  and  yet  preserve  the  quantity  of  information  displayed.  (The  ASV  as  well 
as  any  other  computer  based  navigation  system  could  certainly  use  a  more 
complex  classification  system.  Limiting  the  classification  system  in  favor  of  the 
display  significantly  reduces  the  apparent  range  of  information  which  has  been 
mathematically  derived  about  the  scene.  For  example,  a  line  of  valley  pixels 
oriented  up  and  down  a  slope  could  be  interpreted  as  a  gully  whereas  the  same 
line  oriented  across  a  slope  would  be  a  trail  or  road.  The  differentiation  between 
these  two  terrain  types  is  obviously  important  for  an  autonomous  vehicle.  The 


generation  of  this  kind  of  information  and  its  display  is  a  subject  for  further 
study.) 

In  deciding  on  the  scheme  for  assigning  the  symbols  to  the  classes  of 
terrain,  two  primary  consideration  were  taken  into  account.  The  first  was 
mobility.  As  has  been  stated  previously,  it  is  felt  that  an  autonomous  vehicle 
should  be  primarily  concerned  with  identifying  level  terrain  for  it’s  mobility.  In 
general,  level  terrain  requires  the  least  expenditure  of  energy  to  traverse  and 
therefore  is  more  desirable  in  terms  of  an  energy  cost  function. 

Secondly,  an  autonomous  system  will  not  be  able  to  navigate  adequately  if 
it’s  knowledge  of  it’s  environment  is  strictly  in  terms  of  single  voxels.  The  voxels 
need  to  be  aggregated  into  large  scale  terrain  features  in  order  for  an  autonomous 
vehicle  to  perform  path  planning.  The  knowledge  of  where  the  top  of  a  hill  is  is 
not  as  good  as  knowing  the  lateral  extent  of  that  hill. 

Most  major  terrain  features  seem  to  originate  at  a  level  point,  a  point  with 
zero  or  near  zero  slope.  This  might  be  a  hilltop,  the  bottom  of  a  depression,  the 
floor  of  a  valley,  or  the  top  of  a  ridge  line,  for  example.  Therefore,  an  aggregation 
process  would  probably  use  points  with  essentially  zero  slope  as  seeds  from  which 
to  "grow"  major  terrain  feature.  A  peak  voxel,  with  zero  slope,  which  is 
surrounded  by  safe  or  unsafe  slope  voxels  would  aggregate  into  a  hill. 

In  writting  this  simulation,  major  emphasis  was  given  to  voxels  which  have 
approximately  zero  slope,  that  was  levtl  voxels.  It  was  felt  that  the  maximum 

amount  of  knowledge  possible  about  these  voxels  should  be  displayed  even  if  it  is 
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at  the  expense  of  the  level  of  the  displayed  knowledge  about  the  safe  and  unsafe 
slope  voxels.  Accordingly,  a  symbol  that  can  be  somewhat  associated  with  each  of 
the  secondary  terrain  types,  is  assigned  to  each  voxel  whose  primary  type  is  level. 
All  other  voxels  which  are  classified  as  safe  slope  voxels  are  assigned  a  symbol  of 
"/"  and  unsafe  slope  voxels  are  assigned  the  symbol  "U".  TABLE  8  is  a  table  of 
the  possible  primary  and  secondary  classifiers  and  their  associated  symbols. 


TABLE  8. 

TERRAIN  SYMBOLS  BY 
VOXEL  CLASSIFICATION 


Secondary' 

Classifier 

Level 

Primary  Classifiei 

Safe  Slope 

Unsafe  Slope 

Flat  (Planar) 

- 

/ 

U 

Ridge 

/ 

U 

Saddle 

S 

/ 

U 

Peak 

(up  arrow) 

/ 

U 

Valley 

V 

/ 

U 

Depression 

(down  arrow) 

/ 

U 

Pass 

P 

/ 

U 

Due  to  difficulties  in  producing  the  symbols  for  the  chart,  a  peak  is  presented  as 
an  up  arrow  on  the  symbolic  map,  while  a  depression  is  displayed  as  a  down 
arrow,  and  a  ridge  and  valley  will  be  an  upside  down  and  right  side  up  "v" 
respectively.  In  conclusion,  TABLE  8  can  be  thought  of  as  a  kind  of  legend  for 
reading  the  symbolic  terrain  maps  in  the  next  section. 


E.  SUMMARY 

In  this  chapter,  the  reasons  for  selecting  LISP  as  the  language  in  which  to 
write  this  simulation  were  discussed.  LISP  is  a  rich  and  versatile  language  having 
capabilities  most  other  languages  do  not  have.  One  of  its  primary  strengths 
besides  its  list  processing  functions  is  it’s  ability  to  absorb  foreign  functions.  It 
was  exactly  this  capability  which  allowed  the  program  to  generate  its  results  in 
pictorial  form  instead  of  as  reams  of  numeric  data.  In  this  vein,  a  short  discussion 
was  also  presented  to  acquaint  the  reader  with  the  graphics  capabilities  of  the  ISI 
Optimum  V  Workstations  on  which  the  simulation  was  run. 

The  program  itself  was  also  discussed.  One  of  the  other  strengths  of  LISP 
is  that  the  programmer  is  able  to,  indeed  expected  to,  customize  the  language  to 
his  needs.  The  entire  main  body  of  the  program  consists  of  four  calls  to  one 
function  followed  by  one  call  the  another.  This  is  possible  because  of  the  unique 
structure  of  LISP  wherein  the  programmer  can  create  functions  at  will  and  intern 
them  on  the  oblist  so  that  they  may  be  accessed  just  as  any  other  LISP  function. 

Finally,  the  reasons  for  selecting  the  symbols  with  which  to  display  the 
information  produced  were  presented.  For  this  simulation,  it  was  felt  that  the 
level  terrain  contained  the  most  significant  knowledge  and  so  the  displays  are 
biased  in  that  direction.  The  safe  and  unsafe  slope  terrain  is  displayed,  but  other 
terrain  characteristics  have  to  be  accessed  by  referencing  the  list  of  voxels  named 


"rclass  map"  and  "nclass  map". 


In  the  next  chapter,  the  results  of  the  simulation  will  be  presented.  Each 
artificial  terrain  sample  that  was  created  was  analyzed  a  number  of  times  under 
both  noiseless  and  noisy  conditions.  Also,  for  each  different  scene,  the  two 
thresholds  which  effect  the  results,  the  slope  and  curvature  thresholds,  were 
changed  and  then  the  scene  was  analyzed  again.  The  results  are  both  interesting 
and  surprising  and  show  the  value  of  this  technique  of  terrain  analysis  for 


V.  EXPERIMENTAL  RESULTS 

A.  INTRODUCTION 

In  Chapter  3,  a  mathematical  algorithm  was  developed  to  describe  the 
behavior  of  a  surface  as  modeled  by  a  quadratic  least  squares  fit  to  a  3x3  window 
about  a  point.  The  mathematical  model  was  validated  for  a  3x3  window,  a  case 
which  implies  only  ndividual  point  size  terrain  features.  Except  for  very  large 
scale  images,  most  terrain  features  will  not  be  limited  in  scope  to  individual 
voxels.  For  the  local  navigation  problem,  in  particular,  the  scope  of  an  image  will 
be  such  that  significant  terrain  features  will  span  over  a  number  of  pixels. 
Therefore  it  is  important  to  test  the  algorithm  against  a  number  of  different 
artificial  scenes  to  see  exactly  how  well  it  can  identify  different  major  terrain 
fe..rures.  For  this  test,  the  algorithm  still  fits  a  surface  to  only  a  3x3  window,  but 
it  does  it  to  each  voxel  in  the  scene.  (Along  the  edge,  the  voxels  do  not  have 
enough  neighbors  for  the  algorithm  to  work  so  they  are  immediately  classified  as 
"edge"  voxels.)  Individual  terrain  features  should  stand  out  if  the  method  is  to  be 
of  value  to  an  autonomous  vehicle. 

In  this  chapter,  the  algorithm  is  tested  on  5  different  random  terrain 
scenes.  Each  scene  is  a  square  30  voxels  by  30  voxels  in  size.  These  scenes  range 
from  a  very  benign  scene  with  two  major  features  to  very  complex  scenes  with  a 


large  number  of  individually  identifiable  features.  Each  scene  was  analyzed  six 


different  times,  with  one  or  the  other  of  the  thresholds  which  effect  the 
performance  of  the  algorithm,  curvthreshold  and  slope  limit.  being  changed. 
Lastly,  each  sequence  of  tests  was  performed  against  a  scene  with  no  noise  and 
then  against  the  same  scene  with  a  random  noise  signal  impressed  on  top  of  the 
terrain  signal.  The  random  noise  signal  is  a  random  integer  in  [  -1,  0  .  1  ].  Each 
random  number  was  multiplied  by  0.1  to  produce  a  noise  level  of  xO.l  distance 
units,  where  1  distance  unit  is  the  distance  between  two  voxels  since  the 
algorithm  has  been  non-dimensionalized. 

The  symbolic  results  of  the  5  runs  are  presented  as  Figure  7.  a-n  through 
Figure  11.  a-n.  The  first  two  figures  of  each  run  are  the  noiseless  random  terrain 
and  the  noisy  random  terrain.  The  next  10  figures  are  the  noiseless  and  noisy 
scenes  as  analyzed  by  the  algorithm  with  the  curv  threshold  being  reduced  each 
time.  The  last  two  figures  of  each  run  axe  the  analyzed  noiseless  and  noisy  scenes 
in  which  the  slope  limited  threshold  has  been  increased.  In  all  the  runs  the 
threshold  of  separation  between  safe  slope  and  unsafe  slope  has  been  left  at  30°. 
The  choice  of  30°  was  purely  arbitrary  and  not  based  on  the  operating 
characteristics  of  any  vehicle  in  particular. 

B.  GENERAL  ANALYSIS 

If  the  assumption  is  made  that  an  autonomous  vehicle  utilizing  this 
algorithm  is  capable  of  traversing  over  terrain  which  is  classified  as  either  level  or 
safe  slope,  then  this  algorithm  adequately  analyzes  each  scene  for  the  local 
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Figure  7.  c. 
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Figure  7.  d. 

Curv  threshold  =  0.2 
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Figure  7.  e. 

Curv  threshold  =  0.15 
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Slope  limit  =  2° 
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Figure  7.  f. 

Curv  threshold  =  0.15 
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Figure  7.  g. 
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Figure  7.  h. 

Curv  threshold  =  0.10 


Noisy  Terrain 
Slope  limit  =  2° 
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Figure  7.  i.  Clean  Terrain 

Curv_threshold  =  0.05  Slope  limit  =  2° 


Figure  7.  j.  Noisy  Terrain 

Curv_threshold  =  0.05  Slopelimit  =  2° 
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Figure  7.  m. 
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Figure  7.  n. 

Noisy  Terrain 

Curv  threshold  =  0.05 

Slope  limit  =  10° 
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Figure  8.  c. 
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Figure  8.  d. 
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Figure  8.  e. 

Curv  threshold  =  0.15 
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Slope  limit  =  2 
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Figure  8.  f.  Noisy  Terrain 

Curvthreshold  =  0.15  Slope  limit  =  2 


Figure  8.  g. 

Curv  threshold  =  0.10 


Clean  Terrain 
Slope  limit  =  2C 
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Figure  8.  h.  Noisy  Terrain 

Curv  threshold  =  0.10  Slope  limit  =  2° 
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Figure  8.  j.  Noisy  Terrain 

Curv_threshold  =  0.05  Slope  limit  =  2° 
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Figure  8.  n.  Noisy  Terrain 
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Figure  9.  d. 
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Figure  9.  f. 
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Figure  9.  g. 
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Curv  threshold  =  0.10  Slope  limit  =  2 
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Figure  9.  h. 

Curv  threshold  =  0.10 


Noisy  Terrain 


Slope  limit  =  2 
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Figure  9.  i, 

Curv  threshold  =  0.05 
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Slope  limit  =  2 
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Figure  9.  j. 

Curv  threshold  =  0.05 


Noisy  Terrain 
Slope  limit  =  2 
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Figure  9.  k. 

Curv  threshold  =  0.01 


Clean  Terrain 
Slope  limit  =  2° 
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Figure  9.  1. 

Curv  threshold  =  0.01 


Noisy  Terrain 
Slope  limit  =  2° 
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Figure  9.  m. 

Curv  threshold  =  0.05 


Clean  Terrain 
Slope  limit  =  li 
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Figure  9.  n. 

Curv  threshold  =  0.05 


Noisy  Terrain 
Slope  limit  =  1 
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Figure  10.  c. 
Curv  threshold  =  0.2 


Clean  Terrain 
Slope  limit  =  2 
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Figure  10.  d.  Noisy  Terrain 

Curv  threshold  =  0.2  Slope  limit  =  2 
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Figure  10.  e.  Clean  Terrain 

Curv  threshold  =  0.15  Slope  limit  =  2° 
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Figure  10.  f.  Noisy  Terrain 

Curv  threshold  =  0.15  Slope  limit  =  2° 
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Figure  10.  g. 

Curv  threshold  =  0.10 


Clean  Terrain 
Slope  limit  =  2 
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Figure  10.  h.  Noisy  Terrain 

Curv  threshold  =  0.10  Slope  limit  =  2 
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Figure  10.  i.  Clean  Terrain 

Curv  threshold  =  0.05  Slope  limit  =  2° 
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Figure  10.  j, 

Curv  threshold  =  0.05 


Noisy  Terrain 
Slope  limit  =  2C 
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Figure  10.  k. 

Curv  threshold  =  0.01 


Clean  Terrain 
Slope  limit  =  2° 
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Figure  10.  1. 

Curv  threshold  =  0.01 


Noisy  Terrain 
Slope  limit  =  2C 
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Figure  10.  m. 
Curv  threshold  =  0.05 


Clean  Terrain 
Slope  limit  =  10° 
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Figure  10.  n. 

Curv  threshold  =  0.05 


Noisy  Terrain 
Slope  limit  =  10° 
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Figure  11.  d. 
Curv  threshold  =  0.2 


Noisy  Terrain 


Slope  limit  =  2 
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Figure  11.  e. 

Curv  threshold  =  0.15 


Clean  Terrain 
Slope  limit  =  2° 
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Figure  11.  f. 

Curv  threshold  =  0.15 


Noisy  Terrain 
Slope  limit  =  2° 
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Figure  11.  g. 

Curv  threshold  =  0.10 


Clean  Terrain 
Slope  limit  =  2C 


EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE 

EUUUUUUUUl.'UULKJUU  UUUUUUUE 

EUUUUUUUUU  UUUUU  U'JU  'UUUUUUUE 

EUUUUUUUUUUUUU  U  U  U  U  U  UUU  U  UUE 

E  U  U  U  U  U  EM.'  U  UUUUU  U  U  U  U  U  U  U  U  U  U  E 

EU UUUUUUU  U  UUU  UUUUU  •  UUUUE 

EUUUU  U  U  U  U  U  U  UU  UUUUE 

EUUU.  u u  U UUUUU  U  U  U  E 

£  U  ■UUUUUU  UUUUUUU  U  U  E 


'  U  U  U  U  U  U  U  U  UUUUUUU 
U  U  U  U  U  U  U  U  U  U  U  U  UL'UU  u  u 

U  U  U  U  U  U  U  UUUUUU  UUUUU  u  u 
UUUUUUUUUUUUUU  UUUUUU 

UUUUUUUUUUUUUU  -  UUUUU u  u 
U  U  UUUUU  UUUUUUU  U  UUUU  U  L' 
U  U  U  IJ  UUUUUUUU  U  U  U  UUU  u  u 

U  UUUUUUUUUUUUU U U  '  u u u 

UU UUUUUUUUUUUUUU 
u  u  u  u  u  u  u  u  u  u  u  u  u  u  u  u  u 
U  UUUUUU  UUUUUU  U  U  U  U  L’ 
u  u u u u  u u U UUUU  U UUUU  u u 
U  UUU  U  UUUU  UUUU  L'  u  u 
U  UUUU  U  U  U  U  U  LI  u  u  u 
UUUUUUUUUUUUU 
U  u  U  U  U  U  U  U  L'  U  U 


u  u  u  u 

u  U  U  u  U  L1  U 
U  U  U  U  U  L'  U 
IJ  U  U  U  U  U  L*  I.'  U 


U UUUUUU u u 
UUUUU UUUU 
U  UUUUUU  U 
U  tj  UUUU  U 


Figure  11.  h. 

Curv  threshold  =  0.10 


Noisy  Terrain 
Slope  limit  =  2C 
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Figure  11.  i.  Clean  Terrain 

Curv  threshold  =  0.05  Slope  limit  =  2 
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Figure  11.  j.  Noisy  Terrain 

Curvthreshold  =  0.05  Slope  limit  =  2 
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Figure  11.  1. 
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Figure  11.  m. 

Curv  threshold  =  0.05 
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Figure  11.  n. 

Curv  threshold  =  0.05 


fNoisy  Terrain 
Slope  limit  =  1 


navigation  problem.  In  Figure  7  all  the  terrain  is  either  level  or  safe  slope  and 
traversable.  Figure  8  is  an  interesting  example  in  that  a  second  safe  route,  other 
than  the  route  up  the  pass  between  the  two  peaks,  although  not  obvious  in  the 
terrain  picture,  is  easily  seen  in  the  symbolic  map.  The  second  route  is  along  the 
right  side  of  the  scene  and  then  up  the  side  of  the  right  peak.  The  other  scenes  are 
not  particularly  interesting  except  for  the  fact  that  in  each  case  the  algorithm 
produces  a  symbolic  map  which  can  be  used  to  effectively  pick  out  safe  routes  for 
the  local  navigation  problem  or  which,  in  the  case  of  Figure  II,  will  allow  the 
autonomous  vehicle  to  decide  to  not  traverse  the  scene. 

C.  THE  EFFECT  OF  CURVTHRESHOLD 

The  first  of  the  two  thresholds  to  be  examined  in  this  sequence  of  tests  is 
the  threshold  defined  as  curv threshold  in  the  program.  It  is  the  threshold 
between  a  "significant"  and  an  "unsignificant"  eigenvalue  of  the  Hessian  matrix. 
This  can  be  though  of  most  easily  in  terms  of  the  separation  between  significant 
and  unsignificant  curvature  on  the  surface  fitted  to  the  window.  If  the  curvature 
is  below  the  threshold  value,  it  is  treated  as  a  zero  value  when  trying  to  classify 
terrain  as  per  TABLE  2.  For  all  of  the  runs,  the  value  of  curv  threshold  was 
varied  using  the  following  5  values:  0.2.  0.15.  0.1.  0.05.  and  0.01. 

Run  1,  Figure  7.  a-n.  is  the  most  illustrative  of  the  effects  of  the  change  in 
curv  threshold.  In  the  other  runs  the  slope  limit  over  powers  curv  threshold.  At 
the  largest  value  for  curv  threshold,  in  run  1.  there  is  very  little  detail  available. 


Specifically,  in  the  noiseless  terrain,  all  the  level  terrain  shows  up  as  being  fiat  or 
planar:  i.e.  without  any  curvature.  It  is  not  until  the  curv  threshold  decreases  to 
0.1  that  any  detail  starts  to  surface.  The  small  peak  at  the  upper  left  center  starts 
to  show  as  a  ridge  and  the  low  area  at  the  upper  right  is  classified  as  valley  pixels. 

When  the  curv  threshold  is  decreased  again  to  0.05.  a  significant  amount 
of  detail  about  the  scene  is  brought  out.  One  of  the  pixels  on  the  small  hilltop  is 
classified  as  a  peak  terrain  type.  A  hilltop  is  expected  to  have  peak  types 
surrounded  by  safe  or  unsafe  slopes.  The  flat  areas  are  surrounded  by  valley 
terrain  cells  for  the  most  part.  A  few  of  the  voxels  were  classified  as  ridges.  An 
area  which  starts  to  curve  up  in  one  direction  while  remaining  relatively  flat  in 
the  other  direction  will  have  a  positive  eigenvalue  and  so,  using  TABLE  2.  it 
should  be  expected  that  low  spots  would  have  valley  types  surrounding  them. 

When  the  curv  threshold  is  reduced  to  0.01  many  of  the  valley  and  ridge 
terrain  types  change  to  passes  and  saddles.  The  difference  between  ridges  and 
valleys  and  passes  and  saddles  is  that  the  smaller  eigenvalue  in  magnitude  is  zero 
for  the  ridges  and  valleys  while  the  eigenvalue  is  significant  for  the  pass  and 
saddle.  Both  on  the  dividing  line  around  the  low  areas  and  on  the  peak,  the 
voxel’s  which  have  a  pass  and  saddle  classification  should  probably  be  considered 
in  error. 

Figure  9.  a-n  is  the  only  other  scene  which  has  a  significantly  sized  level 
area  in  which  the  effect  of  curv  threshold  can  be  examined.  This  area  is  sparse 

and  without  detail.  As  the  curv  threshold  is  decreased  the  voxels  take  on  a 
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random  set  of  classifications.  In  an  area  with  little  detail,  a  random  set  of 
classifications  could  be  expected. 

In  all  the  other  scenes,  the  artificial  terrain  has  a  significant  slope  to  it. 
effectively  burying  the  effect  of  curv  threshold.  From  the  two  cases  examined, 
though,  is  seems  that  the  curv  threshold  returns  the  best  information  at  a  value 
between  0.1  and  0.05.  It  may  be  that  upon  a  more  detailed  examination  an 
optimum  value  for  the  curvthreshold  may  be  determinable.  It  seems  more  likely 
that  the  curv  threshold  will  have  to  be  varied  up  and  down  between  about  0.15 
and  0.05.  depending  upon  the  overall  flatness  of  the  scene,  to  generate  good 
knowledge  about  the  terrain. 

D.  THE  EFFECT  OF  SLOPE  LIMIT 

To  investigate  the  effect  of  altering  the  threshold,  slope  limit,  the 
curv  threshold  was  set  to  a  value  of  0.05  which  seemed  to  be  a  reasonable  level 
from  the  previous  section  and  the  slope  limit  was  changed  from  2°  to  10°.  The 
effect  of  this  was  far  more  pronounced  than  the  effect  encountered  in  altering 
curvthreshold. 

In  Figure  8,  not  much  detail  of  the  scene  was  evident  until  the  slopelimit 
was  increased.  When  it  was  increased,  the  pass  between  the  two  peaks  becomes 
easier  to  identify.  The  voxels  between  the  two  hills  are  classified  as  valleys  and 
depressions  for  the  most  part.  A  valley  or  depression  with  unsafe  slope  on  either 
side  of  it  would  be  a  good  way  of  describing  a  pass.  The  rightmost  hill  is  broken 


out  as  a  ridge  voxel.  If  the  curv  threshold  were  reduced  the  voxel  might  show  up 
as  a  peak  cell.  It  is  also  interesting  to  note  that  the  area  on  the  lower  left  side, 
where  the  climbing  slope  breaks  back  down,  is  delineated  with  a  series  of  pass  and 
saddle  voxels.  If  the  curv_threshold  were  increased  in  this  area,  most  of  the 
climbing  slopes  would  probably  show  up  sis  a  ridge  line.  Given  this  information, 
an  autonomous  vehicle  would  probably  identify  the  area  passing  between  the  two 
unsafe  areas  as  some  sort  of  a  pass  and.  therefore,  a  usable  path. 

In  Figure  9.  the  peak  is  finally  broken  out  when  the  slopelimit  is  increased 
to  10°  as  is  the  depression  to  the  right  side  of  the  scene.  The  ridge  line  in  the 
upper  right  side  is  broken  out  with  a  voxel  being  classified  as  a  ridge  type  along 
with  a  series  of  saddle  and  pass  types.  If  the  curv  limit  were  increased  then  these 
cells  would  probably  degrade  to  a  sequence  of  ridge  type  cells. 

Though  more  detail  is  brought  out  in  Figure  10  with  the  increase  in  the 
slope  limit,  the  increase  is  not  large  enough  to  bring  out  significant  detail.  The 
only  small  increase  in  the  detail  is  along  the  valley  in  the  lower  left  corner  in 
which  a  few  of  the  voxels  become  valley  and  pass  type  cells.  An  autonomous 
vehicle  transiting  this  area  would  probably  be  able  to  analyze  the  scene  enough  to 
navigate  it  based  on  the  secondary  characteristics  of  the  safe  slope  cells,  which  on 
this  type  of  symbolic  map  is  not  displayed. 

In  Figure  11,  the  change  in  slope  limit  brings  out  considerable  detail.  The 
depression  in  the  lower  left  corner  shows  up  as  a  depression  voxel  with  a  few  pass 


type  cells  in  the  neighborhood.  The  pass  in  the  upper  right  corner  also  shows  up 


as  a  line  of  valley,  pass  and  depression  voxels.  The  ridge  line  running  from  the 
lower  right  corner  towards  the  upper  left  corner  also  shows  up  well.  Along  this 
line  are  a  number  of  peak  cells  and  saddle  cells  with  a  few  ridge  cells.  With  a  good 
set  of  rules,  an  autonomous  vehicle  could  probably  aggregate  that  sequence  of 
voxels  into  the  ridge  line  that  it  is. 

An  increase  in  the  slope_lim.it  will  not  always  bring  out  an  effective 
increase  in  detail,  though.  Figure  7  is  a  good  example  of  how.  for  a  human 
operator  at  least,  the  resultant  increase  in  displayed  information  with  an  increase 
in  slope  limit  can  actually  increase  the  confusion.  At  a  limit  of  2°  a  nice  amount 
of  detail  is  evident.  But  at  a  limit  of  10°  there  is  too  much  detail  for  a  human 
operator  to  absorb  efficiently,  though  it  should  be  remembered  that  a  computer 
can  effectively  utilize  much  larger  amounts  of  data  than  can  a  human.  It  is 
interesting  to  note,  however,  that  whereas  the  peak  in  the  upper  left  center  area 
had  only  displayed  a  single  peak  voxel  at  a  2°  slope  limit,  there  are  a  number  of 
peak  cells  in  the  area  when  the  slope  limit  is  increased  to  10°.  This  would  be  a 
very  strong  indication  of  an  actual  peak  instead  of  a  false  peak  due  to  noise. 

E.  THE  EFFECT  OF  NOISE 

Each  scene  is  analyzed  under  both  noiseless  and  noisy  conditions.  The 
results  of  both  analyses  are  presented  on  the  same  page  for  easy  comparison.  As  a 
general  observation,  there  are  two  different  levels  at  which  the  noise  effects  the 


analysis.  The  first  level  is  for  the  primary  classification  and  the  second  is  for  the 


secondary  classification. 

1.  The  Effect  On  Primary  Classification 

At  the  level  of  generating  a  primary  classification,  the  effect  of  noise 
is  strongest,  as  would  be  expected,  at  the  boundaries  between  level,  safe-slope  and 
unsafe-slope  areas.  However,  at  a  slope  limit  of  2°  the  level  areas  are  very 
sensitive  to  noise  as  shown  in  Figures  7  d.  and  9  d.  where  the  boundary  as  well  as 
the  interior  of  the  area  is  effected  by  the  noise.  At  a  slope  limit  of  10°.  though, 
the  effect  of  noise  is  mostly  on  the  edges  of  the  area  as  shown  in  Figure  9  n.  In  all 
the  other  runs,  at  both  the  2°  and  10°  slope  limit,  the  noise  has  little  effect  on 
the  safe-slope  and  unsafe-slope  areas  except,  again,  at  the  boundary.  In  general, 
the  boundaries  are  usual'y  displaced  only  by  a  single  voxel.  For  the  purposes  of 
this  algorithm,  this  level  of  noise  is  acceptable  since  the  separation  between  safe- 
slope  and  unsafe-slope  can  be  chosen  conservatively  and  rules  can  be  entered  into 
the  navigation  package  allowing  an  autonomous  vehicle  to  traverse  no  closer  than 
2  (or  whatever  number  is  desired)  voxels  from  a  particular  boundary. 

2.  The  Effect  On  Secondary  Classification 

The  noise  has  the  largest  effect  on  the  secondary  classification  of  a 
voxel.  While,  in  this  algorithm,  the  secondary  classification  of  a  voxel  is  only 
displayed  when  the  voxel  has  a  primary  classification  of  "level",  the  effect  of  noise 
can  still  be  seen  clearly  in  Run  1  and  in  the  final  examples  of  Runs  2  and  5.  In 

particular,  the  noise  alters  the  eigenvalues  enough  so  that  the  algorithm  falsely 
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classifies  a  large  number  of  the  voxels.  In  Run  1.  the  difference  between  the  noisy 
and  noiseless  symbolic  maps  is  significant  with  the  noisy  map  having  a  large 
number  of  random  terrain  types  while  the  noiseless  map  has  few.  if  any.  voxels 
identified  as  other  than  "flat"  (or  planar).  In  Runs  2  and  5  the  noise  changes  the 
classification  of  a  number  of  voxels  but  the  effect  is  not  as  significant  as  when  the 
amplitude  of  the  terrain  is  relatively  flat,  as  is  the  case  for  Run  1.  Considering 
that  a  pass  is  a  subset  of  a  valley  and  a  saddle  is  a  subset  of  a  ridge,  when  looking 
at  the  final  examples  of  Runs  2  and  5.  it  appears  that  the  general  classification  of 
a  voxel  is  changed  little  by  the  noise.  A  number  of  valleys  turn  into  passes  and  a 
number  of  ridges  turn  into  saddles,  but  only  in  a  few  cases  do  valleys  turn  into 
their  direct  opposites,  ridges,  and  peaks  into  their  opposites,  pits.  It  appears  that 
by  averaging  the  terrain  types  in  the  local  area,  the  effect  of  the  noise  can  be 
partially  countered. 

F.  CONCLUSIONS 

In  this  chapter,  the  results  obtained  from  applying  the  algorithm  to  a  large 
area,  instead  of  just  a  single  3x3  window,  were  examined.  The  performance  of  the 
algorithm  was  validated  in  a  previous  chapter  against  point  sized  terrain  features, 
but  few  major  terrain  features  are  point  sized. 

From  the  results  examined  in  this  chapter,  it  appears  that  the  first 
derivative  information  is  of  sufficient  quality  for  an  autonomous  vehicle  to 
navigate.  Level  and  safe  slope  routes  through  the  scenes  were  clearly  displayed 
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and  in  one  case  a  non-obvious  route  was  found  by  the  algorithm.  The  second 
derivative  information  added  significant  detail  to  the  symbolic  map.  Individual 
peaks  and  pits,  valleys  and  ridges,  and  saddles  and  passes  were  highlighted.  This 
information  would  allow  an  autonomous  vehicle  to  not  only  navigate  based  on 
safe  and  unsafe  areas  but  also  based  on  major  terrain  features  such  as  passes 
through  hills,  peak  and  pits,  etc. 

The  performance  of  the  simulation  is  very  dependent  on  three  variables. 
The  variable  first  is  the  slope  limit  which  is  the  cut  off  between  level  and 
safeslope  classified  terrain.  (Another  threshold  exists  between  unsafe  slope  and 
safeslope  which  was  not  examined.)  The  effect  of  this  threshold  was  examined  by 
varying  it  between  2°  and  10°.  At  the  lower  threshold,  a  good  amount  of  detail  in 
the  low  amplitude  run.  Run  1.  was  displayed  but  little  detail  in  the  other  runs 
was  displayed.  At  the  upper  threshold,  the  detail  overpowered  the  first  derivative 
information  in  Run  1  but  was  good  for  the  other  runs.  The  second  variable  is 
curvthreshold  which  separates  the  significant  and  nonsignificant  eigenvalues. 
Generally,  the  lower  the  threshold  the  better  the  information  displayed.  This  was 
true  until  the  threshold  was  lowered  down  to  0.01  at  which  point  the  difference 
between  a  valley  and  a  pass  and  a  ridge  and  a  saddle  became  indistinct.  Also, 
when  the  terrain  was  essentially  flat,  the  lower  threshold  seemed  to  force  another 
classification  on  what  were  really  flat  or  planar  voxels.  Finally,  the  effect  of  noise 
on  the  algorithm  was  examined.  Noise  appeared  to  have  a  smaller  effect  on  the 
terrain  scenes  with  larger  vertical  deviations  while  the  scene  with  little  vertical 
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development.  Run  1,  was  strongly  effected  by  the  small  amount  of  noise 
introduced  into  the  signal.  In  other  words,  the  large  amplitude  terrain  signals 
overpowered  the  noise  while  the  low  power  terrain  signals  were  over  powered  by 
the  noise. 

As  can  be  seen  from  the  results  presented  in  this  chapter,  the  simulation  is 
very  dependent  on  the  two  parameters,  slope_limit  and  curv ^threshold.  Before 
this  algorithm  can  be  effectively  utilized  in  an  autonomous  vehicle,  appropriate 
values  for  these  thresholds  will  have  to  be  determined.  In  trying  to  select  values 
for  the  thresholds,  it  may  be  discovered  that  a  form  of  adaptive  thresholding 
based  on  the  power  in  the  terrain  signal  will  actually  be  most  effective.  How  to 
handle  the  noise  is  a  more  significant  problem,  though,  requiring  more  study.  In 
general  terms,  the  noise  behaves  as  noise  does  in  an  electrical  system, 
overpowering  low  power  signals  and  being  overpowered  by  large  amplitude 
signals.  If  normal  low  pass  filtering  concepts  are  employed  to  help  control  the 
noise  (such  as  averaging  the  voxel  types  in  a  local  area),  then  individual  hazards 
(such  as  a  rock  in  a  vehicle’s  path  which  would  show  up  as  a  single  peak  voxel) 
may  be  filtered  out.  Clearly,  noise  is  a  difficult  problem  requiring  careful  analysis 
and  a  much  more  work  than  the  scope  of  this  thesis  permitted. 
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VI.  SUMMARY  AND  CONCLUSIONS 


Any  autonomous  vehicle,  whether  it  is  a  land  vehicle  walking  the  hills  of 
Colorado  or  a  space  probe  exploring  the  planet  Mars,  requires  the  ability  to  sense 
it's  environment.  If  the  vehicle  operates  at  a  distance  from  it's  target,  such  as  the 
Voyager  probes  sensing  Saturn  at  thousands  of  kilometers  or  a  Cruise  Missile 
flying  at  a  predetermined  altitude,  then  the  sensors  do  not  require  extremely  high 
resolution.  Sensors  operating  in  the  IR  region,  Radar,  and  Ultrasound  may  be 
good  enough  in  those  cases.  However,  if  the  vehicle  operates  in  very  close 
proximity  to  a  hostile  environment,  such  as  by  walking  on  it.  then  it  is 
advantageous  to  sense  the  environment  with  high  resolution.  This  implies  utilizing 
sensors  operating  in  the  optical  regime. 

To  date,  most  research  in  optical  wavelength  sensing  has  concentrated  on 
two-dimensional  processing  to  generate  intrinsic  images  or  object  edges.  Little 
research  had  been  conducted  on  three-dimensional  optical  sensing  for  object 
classification,  which  is  more  useful  when  considering  the  autonomous  vehicle 
navigation  problem,  until  the  work  done  by  Olivier  at  Ohio  State  University.  The 
purpose  of  this  thesis,  has  been  to  do  a  detailed  mathematical  development  of  the 
formulas  utilized  by  Olivier  and  then  to  write  a  simulation,  based  on  those 
formulas,  which  attempts  to  analyze  large  random  geographic  terrain  scenes  as  a 
proof  of  the  concept  of  terrain  classification  for  autonomous  vehicle  navigation. 


A. 


RESEARCH  CONTRIBUTIONS 


In  this  thesis,  a  number  of  advances  have  been  made  over  the  work 
originally  publisher  by  Olivier  in  his  masters  thesis.  The  primary  advance  is  the 
mathematical  derivation  of  two  new  terrain  types.  These  are  the  pass  and  the 
saddle.  The  pass  is  a  subtype  of  the  valley  in  which,  instead  of  the  smaller 
eigenvalue  of  the  Hessian  matrix  being  zero,  it  is  negative  while  the  larger 
eigenvalue  is  positive.  The  saddle  is  a  subset  of  the  ridge  in  which,  instead  of  the 
smaller  eigenvalue  being  zero,  the  smaller  eigenvalue  is  positive  while  the  larger 
eigenvalue  is  negative.  These  two  new  terrain  classifications  resulted  from  a 
mathematical  study  of  the  theory  of  quadratic  surfaces.  When  combined  with  the 
five  other  types  previously  defined  in  Olivier's  work,  a  complete  set  of  all  the 
possible  terrain  classifications,  which  can  be  made  based  solely  on  the  values  of 
the  eigenvalues  of  the  Hessian  matrix,  is  formed. 

Secondly,  a  complete  mathematical  derivation  of  the  formulas  on  which 
the  classification  of  terrain  voxels  is  based  was  presented.  The  formulas  are 
dependent  on  the  sampling  points  being  regularly  spaced  about  the  voxel  of 
interest  at  =;  one  distance  unit.  The  particular  distance  unit  is  not  important;  it 
can  be  meters,  centimeters,  or  kilometers  depending  on  the  scale  of  the  scene 
being  analyzed.  The  effect  of  the  regular  spacing  of  the  sampling  points  is  the 
elimination  of  a  number  of  terms  in  the  six  equations  resulting  from  a  least 
squares  error  analysis  of  a  quadratic  equation  describing  a  three  dimensional 
surface.  When  the  six  equations  are  simplified,  through  the  elimination  of  those 


terms,  five  equations  are  derived  which  allow  the  constants  of  the  quadratic  to  be 
found.  These  constants  can  then  be  used  to  easily  calculate  the  gradient  and 
eigenvalues  of  the  Hessian  matrix  for  the  least  square  surface,  which  in  turn,  are 
the  values  which  are  the  basis  of  the  classification  scheme. 

Finally,  the  algorithm,  including  the  two  new  terrain  classification  types, 
was  implemented,  in  LISP  (one  of  the  two  main  languages  used  today  in  artificial 
intelligence  research),  and  ran  against  terrain  scenes  with  relatively  realistic  major 
terrain  features.  The  output  of  the  simulation  was  designed  to  graphically  display 
the  capability  of  the  algorithm  to  produce  symbolic  results  which  could  be  used 
by  an  autonomous  vehicle  to  find  safe  paths  through  the  scene  as  well  as  identify 
major  terrain  features  for  the  purposes  of  navigation  planning.  From  the  results  of 
the  five  simulation  runs,  it  is  clear  that  the  algorithm  is  capable  of  analyzing  large 
scale  scenes.  It  is  also  clear  that  the  performance  of  the  algorithm  is  very 
dependent  on  the  proper  selection  of  appropriate  thresholds  for  the  roughness  of 
the  terrain  being  analyzed.  If  the  slope_limit  is  too  low  or  the  curv_limit  is  too 
high,  the  amount  of  knowledge  generated  by  the  simulation  is  degraded.  If  the 
two  thresholds  are  set  to  be  too  sensitive,  the  information  generated  is  of 
questionable  value.  The  effect  of  noise  on  the  performance  of  the  simulation  was 
also  examined  and  found  to  be  significant  depending  on  the  "power"  in  the 
terrain  signal. 

Another  major  result  of  this  work,  though  not  related  specifically  to  the 
algorithm,  is  the  generation  of  a  number  of  routines  to  handle  graphics  in  LISP 


on  the  1SI  Optimum  Y  Workstation.  Franz  LISP,  the  variant  of  LISP 
implemented  on  the  workstations,  has  no  inherent  graphics  capabilities.  This  is  a 
major  limitation  considering  the  outstanding  graphics  capabilities  of  the 
workstations.  By  "importing"  these  graphics  routines  into  LISP  (through  the 
LISP  cfasl  command),  new  LISP  functions  can  be  generated  which  will  allow 
Franz  LISP  to  utilize  most  of  the  graphics  capabilities  of  the  workstations. 

B.  RESEARCH  EXTENSIONS 

This  simulation  was  developed  primarily  as  a  proof  of  concept  of  the 
usefulness  of  the  idea  of  terrain  classification  based  on  range  images.  As  such  it 
has  been  very  successful.  There  are  a  number  of  problems,  though,  with  the 
simulation  which  require  further  study  before  the  algorithm  can  be  employed  in 
an  autonomous  vehicle.  First,  the  algorithm  is  relatively  slow.  This  simulation  is 
installed  as  an  interperted  set  of  routines  and  has  not  been  compiled  successfully. 
As  such,  it  takes  about  20  minutes  to  complete  one  scene  analysis.  This  is 
prohibitively  long  for  an  autonomous  vehicle.  Therefore,  the  classification  portion 
of  the  simulation  needs  to  be  compiled  on  a  LISP  machine  or  possibly  installed  in 
machine  language  before  it  can  be  used  in  real  time.  It  may  also  be  possible  to 
significantly  improve  the  speed  of  the  simulation  by  writing  the  mathematically 
intensive  portions  of  the  algorithm  in  Fortran  or  Pascal,  compiling  them 
separately,  and  then  accessing  them,  in  the  LISP  environment,  through  the  cfasl 


command. 


.  1 


Another  important  area  for  further  study  is  attempting  to  utilize  the  voxel 
classifications  in  aggregating  individual  voxels  into  major  terrain  features.  It  is 
easy  enough  for  a  human  to  identify  ridge  lines  as  a  series  of  ridge  voxels  and  hills 
as  peak  voxels  surrounded  by  safe-slope  voxels,  but  this  a  capability  the  algorithm 
does  not  now  possess.  By  examining  the  symbolic  results  in  Chapter  5.  it  is 
obvious  that  aggregation  will  involve  more  that  just  combining  like  terrain  types. 
Depending  on  the  thresholding,  for  example,  valley  voxels  can  be  displayed  as 
either  valleys,  passes,  or  even  flat  or  planar  voxels.  The  same  is  true  of  all  the 
other  terrain  types.  Also,  there  are  a  number  of  other  factors  that  need  to  be 
considered  when  attempting  to  aggregate  and  identify  major  terrain  features.  An 
example  is  the  angle  between  the  principal  eigenvector  and  the  gradient  of  a 
voxel.  As  was  briefly  discussed  before,  a  valley  line  oriented  up  a  slope  is  probably 
just  a  ravine,  whereas,  if  it  is  oriented  perpendicular  to  the  slope  it  could  be  a 
path.  The  possible  variation  between  slope,  the  eigenvalues,  the  eigenvectors,  as 
well  as  the  classification  of  nearest  neighbors  leads  to  a  large  number  of 
possibilities.  An  aggregation  scheme  based  on  this  classification  algorithm  is  a 
good  example  of  a  system  which  would  respond  well  to  the  techniques  of  artificial 
intelligence  wherein  a  number  of  general  rules,  based  on  a  some  of  the  factors, 
could  be  devised  to  classify  each  voxel. 

Another  good  way  of  trying  to  aggregate,  instead  of  attempting  to  write 
long  lists  of  rules  which  cover  each  individual  case,  would  be  to  let  a  robot  learn 


what  is  important  when  aggregating  voxels  by  exploring  a  wide  variety  of  terrain 


scenes  and  then  have  the  robot  write  its  own  set  of  rules.  An  actual  robot  in  a 
natural  environment  would  not  be  specifically  required.  Instead  a  computer 
simulation  of  the  dynamics  of  a  robot  or  autonomous  vehicle  could  be  written 
which  would  then  drive  or  walk  across  a  number  of  artificially  created  terrain 
scenes.  By  calculating  the  energy  output  of  the  dynamic  system  ai  d  also 
recording  the  events  which  lead  to  catastrophic  situations  and  then  comparing 
them  to  terrain  classifications  derived  through  this  model,  the  robot  could  "learn" 
to  effectively  utilize  the  information  generated  in  this  simulation.  A  well  written 
simulation  in  LISP  could  potentially  even  write  its  own  set  of  rules  for  both 
aggregation  and  path  planning  which  the  human  could  access  to  see  what  the 
robot  "thinks"  is  important. 

The  effect  of  noise  is  well  worth  further  research.  It  was  only  touched  upon 
in  this  thesis  but  it  is  obvious  that  a  form  of  filtering  must  be  developed  if  the 
algorithm  is  to  be  installed  in  truly  autonomous  vehicles.  A  form  of  low  pass 
filtering,  in  which  neighbors  are  averaged,  may  be  useful  but  it  appears  that  such 
filtering  could  destroy  knowledge  about  small  hazards.  Another  possibility  is  a 
form  of  Kalman  filtering  or  time  averaging  where  random  noise  will  show  up  as 
random  classifications  in  subsequent  sweeps  of  the  sensor,  unlike  a  true  hazard 
which  will  remain  relatively  constant  over  time. 

There  are  numerous  of  other  possible  research  extensions  which  can  evolve 
from  this  simulation.  Adaptive  thresholding  for  the  slopejimit  and 
curv  threshold  and  extending  this  algorithm  to  third  order  or  higher  surfaces  are 


but  a  few.  In  actuality  this  algorithm  has  been  developed  more  as  a  tool  which 


will  hopefully  be  expanded  upon  than  as  an  actual  algorithm  to  be  implemented 
in  a  working  system.  As  such  it  exhibits  a  good  potential  to  be  the  basis  for  giving 
an  autonomous  vehicle  a  reasonable  vision  capability  for  terrain  recognition. 


APPENDIX  A  -  MATRIX  DERIVATION  OF  COEFFICIENTS 


This  is  a  second  derivation  to  find  the  formulas  for  the  coefficients  of  a 
quadratic  least  squares  fit  to  a  three  dimensional  surface.  This  derivation  is  based 
on  matrix  calculus.  First  it  is  necessary  to  assume,  as  in  the  previous  derivation, 
that  the  data  will  be  in  a  regularly  spaced  3x3  window  as  in  TABLE  1.  The 
ordering  of  the  elevation  data  is  as  in  TABLE  9. 


TABLE  9 

ORDERING  OF 
ELEVATION  DATA 
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The  estimating  equation  for  the  elevation  data  is: 

2,  =&,-!-  V.  *  kzV*  *  V,2  +  M.A  *  knV ,2' 

The  error  between  the  real  elevation  z  and  z,  t.  is: 

=  (z  ~  zi) 


(A.l) 


(A.2) 


and  the  sum  of  the  squared  errors, is: 


J 


In  vector  matrix  form.  Equation  A.l  is: 


c  =j  u 
~9il  *^9i6  ^611 

where  the  subscripts  are  the  size  of  the  matrices, 

k  =  (  k2  k2  k4  kS  k6  f 

and 

2  ° 

A  =  (  1  x  y  x  x  u  y~  ). 

t  '  t  I  1*1  t 

The  squared  errors  are  then,  in  matrix  form: 

<Sf  =  (z  -  Z)T  (z  -  z) 

=  (;  -  Ak)T  (z  -  Ak ) 

=  {zT  -  kT  AT)  (z  -  A*) 

=  zT z  -  kT AT z  -  zT Ak  kT AT Ak. 

Letting  the  product: 

^  T 

A  z  =  w 

the  two  middle  terms  can  be  replaced  by 


respectively.  Equations  A.  12  and  A.  13  are  equal  scalar  dot  product  expressions 


and  so  Equation  A.  10  may  be  written  as: 


2kTATz  -  kT A  T Ak. 


(A. 14) 


The  gradient  of  the  error  is: 


= 


dk 

t 


=  -2  ATz  +  2  ATAk  =  0 


(A. 14) 


or 


A  T Ak  =  A  Tz. 


Thus: 


ata 


(A- 15) 


(A. 16) 


or 


=  Bz 


(A.17) 


where 


T 

A  A 


and 
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(A.18) 


1-1  1  1-1  1 


10  10  0  1 

!  i  i  i  i  i  i 

i 

1-1  0  1  0  0 

A  =  1  0  0  0  0  0  (A. 19) 

110  10  0 

1-1-1  1  1  1 

10-1001 
1  1-1  1-1  1 

To  find  the  coefficients  of  the  quadratic  it  is  only  necessary  then  to  calculate  the 
B  matrix  and  then  multiply  it  by  the  z  matrix.  To  calculate  the  B  matrix,  the 
transpose  of  the  A  matrix  must  be  produced  so  that  it  may  be  multiplied  by  the 
A  matrix.  The  transpose  of  the  A  matrix  is: 

111111111 
-10  1-10  1-10  1 

1  1  1  000-1-1-1 

(A. 20) 

101101101 
-10100010-1 
1  1  1  0  0  0  1  1  1 


The  product  of  A  transpose  and  A  is: 


T 

A  A 


(A. 21) 


9  0  0  6  0  6 
!  0  6  0  0  0  0 
0  0  6  0  0  0 
6  0  0  6  0  4 
0  0  0  0  4  0 
6  0  0  4  0  6 

This  product  is  inverted  and  then  multiplied  by  A  transpose.  Equation  A. 20.  The 
inverted  matrix  is: 


UT*)~l 


5/9  0 

0  1/6 
0 

-1  3  0 

0  0 
-1,3  0 


0  1/6 


0  -1/3  0  -1/3 

0  0  0  0 

0  0  0 

0  1/2  0  0 

0  01/4  0 


(A. 22) 


0 


0 


0  1/2 


The  product  of  the  inverted  term  and  A  transpose  is  the  B  matrix,  which  is: 

-1/9  2,9  -1/9  2/9  5/9  2  9  —1/9  2/9  -1/9 

-1/6  0  1/6  -1/6  0  1  6  —1/6  0  1/6 

1/6  1  6  1/6  0  0  0  -1/6  -1/6  -1/6 

1/6  -13  16  16  -1/3  16  16  -  1/3  1/6 

-1  4  0  1/4  0  0  0  1  4  0  —1/4 

1'6  1/6  16  -13  -1/3  -1,3  1/6  1/6  16 


B  = 


(A. 23) 


To  calculate  the  constants  of  the  quadratic  least  squares  surface,  the  above 
matrix,  B.  is  multiplied  by  the  1x9  column  vector  of  elevation  data,  z.  The 
resulting  1x6  column  vector  are  the  6  constants,  through  fcg. 


APPENDIX  B  -  GRAPHICS  FUNCTIONS 


K*X***x**«*9*:¥*xxx:KX**xxx***jtx**jlcxxj|cg( 


TOFILE 

*  / 

This  is  a  C  routine  compiled  separately 
which  is  designed  to  read  a  bitmap  from 
the  ISI  Graphics  machines.  It  is  intended 
to  be  compiled  in  line  in  LISP  environment 
to  provide  graphics  calls.  The  function  * 
is  called  "tofile"  from  LISP.  It  is  called  * 
as  (tofile  bmd  filename)  and  it  will  * 

return  a  t  when  complete.  The  parameters 
passed  to  this  "C"  functions  are:  *, 

*  / 

bmd  =  bitmap  descriptor  number 

*  / 

filename  =  the  file  name  that  the  *' 
bitmap  is  to  be  save  *  ' 

under  * , 


The  way  this  function  works  is  by 
calling  a  C  function  which  reads  the 
bits  in  the  bitmap,  4  at  a  time.  If  the 
first  bit  is  a  1.  then  8  is  added  to  * 

the  counter  bitval.  If  the  next  bit  is  * 
set  then  a  4  is  added  to  bitval.  Etc... 
until  the  value  of  all  4  bits  is  added  * 
to  bitval.  It  then  goes  down  a  list.  * 

comparing  the  value  of  bitval  until 
it  matches  it.  The  letter  or  decimal 
value  that  that  equates  to  in  hexa¬ 
decimal  is  then  printed  to  the  file  * 

under  filename.  The  junk  it  prints  to 
filename  at  the  start  and  finish  of  the 
program  are  control  codes  to  the  laser 
printer  which  allow  the  bitmap  code 
to  be  printed.  It  should  be  noted 
that  the  area  of  the  the  bit  map  that 
is  read  is  576  by  448  bits.  This  is 
one  entire  window  in  this  simulation. 

If  it  is  desired  to  use  this  program  * 

to  read  a  bit  map  of  a  different  size, 
the  area  read  must  be  in  blocks  of  64 
bits  long.  Also,  the  first  numbers  printed 
in  the  file  after  the  control  codes  must 
a  4  digit  decimal  number  which  is  the 
number  of  bits  long  each  row  to  be 


’  printed  is.  After  that,  all  the  other 

*  numbers  are  the  hex  code  of  the  bitmap 
-  *  * 

.  *  * 

/«**«» ***k«k***x****«**xss***z*xx****x***x*k*x*x*x* 


/* 

*  . 

/* 

Dennis  Poulos 

* 

I9 

x  3403 

*  j 

t* 

/* 

26  May  1986 

* 

*  ' 

/*******************»«*******«**********»******»***/ 
★include  <stdio.h> 

★include  <vt.h> 

★  include  <bitmap.h> 

★  include  <math.h> 

FILE  *fopen(),*fp: 


tofile  (bmp, filename) 
char  ‘filename;;: 
int  *bmp; 

{ 

int  temp; 

struct  BMD  *bmd: 

short  iij,k,bitv: 

int  bitval: 

short  boxw,boxh, 

temp  =  *bmp; 

bmd  =  temp: 

fp=fopen(filename."w"); 

fprintfffp."  “PY  *-0F  'GO): 

fprintfffp."  'IJ01500  TT01000); 

fprintfffp," *GOSM11000IGV'PW050IGE"); 

fprintfffp,"  “IP02020); 

boxw=576; 

boxh=448: 

fprintfffp,"  *P"); 

fprintfffp,  "0"); 

fprintff  fp,"%d"  .boxw ) ; 

forfii=0:ii<boxh:ii-r^) 

{ 

forf  j =0y  <  boxw  j  *  =  4 ) 

{bitval  =  bitv  =  0: 

iff(bitv=BM_ReadBit(bmdj,ii))  ==0) 
bitval  ->-=  8; 

iff(bitv=BM_ReadBit(bmd,(j-t-I),ii))  =  =0) 
bitval  —  =  4; 

iff  ( bitv  =  BM_ReadBit(  bmd,  (j  +  2),ii))  ==0) 
bitval  -*-=  2; 

if((bitv=BM_ReadBit(bmd,(j  +  3),ii))  =  =0) 
bitval  ■*  =  1; 


if(  bit  val  =  =  15) 
fprintf(fp."F"): 
else  if(bitval  ==  14) 
fprintf(fp."E"): 
else  if(bitval  ==  13) 
fprintf(fp."D"): 
else  iffbitval  ==  12) 
fprintf(fp."C"): 
else  if(bitval  ==11) 
fprintf(fp."B"); 
else  if(bitval  ==  10) 
fprintf(fp."A"): 

else 

fprintf(fp,"%d".bitval); 

} 

fprintf(fp."0): 

} 

fprintf(fp."  ‘.0G0): 
fprintfjfp,"  *0  'GO-  “PN  ‘-0); 
fclose(fp); 
return: 

} 

/*****9***x*m******** *«******x*«**x*x***x*x***x*»*x  , 

/*  *  ’ 

*  compile  as  * 

'*  cc  -c  tofile.c  -ltools  -lbm  -lvt  * 

»  * 

rx«************«************* ******* ****«*****«****/ 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


OPEN-WINDOW 

This  is  a  fortran  subroutine  compiled 
separately  to  open  a  window  on  the  ISI  * 
Graphics  Machines.  It  is  intended  to  be 
compiled  inline  in  a  LISP  environment  to  * 
provide  LISP  graphics  calls.  The  function  * 
is  called  "opw".  From  LISP  it  is  called  * 
as  (opw  wx  wy  ww  wh  "title")  and  it  will  ’ 
return  the  fdnum  (file  descriptor  number)  * 
when  complete.  The  parameters  passed  to 

the  fortran  function  from  LISP  are:  * 

* 

wx  =  x  coordinate  of  upper  left  corner  * 

* 

wy  =  y  coordinate  of  upper  left  corner  * 

* 

ww  -  window  width  * 

* 

wh  =  window  height  * 

* 

title  =  window  title  (max  31  characters)  * 


*  all  coordinates  are  in  numbers  of  pixels  * 

*  * 


*  * 

*  Dennis  Poulos  * 

*  x  3403  * 

*  15  April  1986  (version  2)  * 

*  * 


integer  function  opw(wx,wy,ww.wh. title) 
external  Open  Window 
integer  Open  Window,  wx.wy,ww,  wh, opw 
character*31  title 

opw  =  Open  Window(wx.wy.ww,wh, title) 

return 

end 

************************************************** 
*  * 

*  compile  as  * 

*  177  -c  openw.f  -ltoolsf  -lbmf  -lvtf  * 

*  * 

************************************************** 
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LINE-DISCIPLINE 

This  is  a  fortran  subroutine  compiled 
separately  to  set  the  line  discipline 
of  the  window  associated  with  the  file 
descriptor  fd.  The  default  line  discipline  * 
of  a  window  is  NTTYDISC.  To  access  the  * 
graphics  capabilities  of  the  ISI  the 
line  discipline  must  be  set  to  TWSDISC. 

The  function  is  called  "slined".  From  LISP  * 
it  is  called  as  (slined  fd)  and  it  will  •* 
return  a  "t"  when  complete.  The  parameters  * 
passed  to  the  fortran  function  from  LISP  * 


*  fd  =  file  descriptor  of  the  desired  * 

*  window  * 

*  * 

*X*XXxxX***X****X*X*x**tt*«X*x****«****x***** ****** 
*  * 

*  Dennis  Poulos  * 

*  x  3403 

*  15  April  1986 

*  * 

integer  function  slined(fd.mode) 
external  SetLineDisc 
integer  fd. slined, mode, SetLineDisc 
slined  =  SetLineDisc(fd.mode) 
return 
end 

*  * 

*  compile  as  * 

*  f77  -c  setldisc.f  -ltoolsf  -lbmf  -lvtf  * 

*  x 
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if 


tx*xxxxxxxx«x*xvx*xxxxxxXxxxxxxxxxxxxx*vxsx*«x*a 


BIT  MAP  ALLOCATE 


This  is  a  foreran  subroutine  compiled 
separately  to  allocate  a  portion  of  memory 
for  the  bit  map.  A  bit  map  is  maintained 
in  memory  and  can  be  displayed  in  a 
window  as  desired.  This  routine  is  intended 
to  be  compiled  inline  in  a  LISP  envir¬ 
onment  to  provide  LISP  graphics  calls. 

The  function  is  called  "bitmapa".  From 
LISP  it  is  called  as  (bitmapa  width  height) 
and  it  will  return  the  bmd  (bit  map  * 

descriptor)  when  complete.  The  parameters 
passed  to  the  fortran  function  from  LISP 


width  =  bit  map  width 


height  -  bit  map  height 


all  coordinates  are  in  numbers  of  pixels  * 


Dennis  Poulos 
x  3403 

18  April  1986 


integer  function  bitmapa(width, height) 
external  bmAllocate 

integer  bmAllocate, bitmapa. width, height 
bitmapa^bm  Allocate)  width,  height) 


compile  as  * 

177  -c  bitmal.f  -Itoolsf  -lbmf  -Ivtf  * 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


ADDRESSING  MODE 

This  is  a  fortran  subroutine  compiled 

separately  to  set  the  display  mode  in  a 

window.  The  two  modes  available  are 

"VT  RELATIVE”  AND  "VT  ABSOLUTE”.  Relative  * 

addressing  means  pixels  are  specified 

relative  to  the  current  pointer  position  * 

while  absolute  addressing  specifies  pixels  * 

relative  to  the  upper  left  hand  corner 

of  the  bitmap.  This  routine  is  intended  * 

to  be  compiled  inline  in  a  LISP  environment  * 

to  provide  LISP  graphics  calls.  The  function  * 

is  called  as  (bmadres  bmd  mode)  and  it  will  * 

return  a  ”t”  when  complete.  The  parameters  * 

passed  to  the  fortran  subroutine  from  * 

LISP  are: 

* 

bmd  =  bitmap  descriptor  number 

* 

mode  =  addressing  mode 

* 

************************************************* 


*  * 

v  Dennis  Poulos  * 

*  x  3403 

*  20  April  1986 

*  * 

************************************************** 

subroutine  bmadres(bmd.mode) 
integer  bmd. mode 
call  bmSet Address(bmd.mode) 
return 
end 


c 


compile  as 

(77  -c  bmadres. f  -Itoolsf  -Ibmf  -Ivtf 


c 


z*****z*zz*z«»z*zz«***» 


c 

c 

c 

c 

c 


SET  BACKGROUND  COLOR 

This  is  a  fortran  subroutine  compiled 
separately  to  set  the  background  color 
of  a  bitmap.  The  background  color  can  * 
be  set  to  either  VT_White  or  VT_BLACK. 

This  routine  is  intended  to  be  compiled  * 
inline  in  a  LISP  environment  to  provide 
LISP  graphics  calls.  The  function  is  called  * 
as  (setbacc  bmd  color)  and  it  will  * 

return  a  "t"  when  complete.  The  * 

parameters  passed  to  the  fortran  subroutine  * 

from  LISP  are:  * 

* 

bmd  =  bitmap  descriptor  number  * 

color  =  VT  White 
or 

VT  Black 

*«***« ******************************************* 


Dennis  Poulos  * 

x  3403 

22  April  1986 

* 


subroutine  setbacc(bmd. color) 
integer  bmd. color 
call  bmSetBColor(bmd. color) 
return 
end 

t***********?************************************* 


c 


compile  as 

f77  -c  setbacc. f  -Itoolsf  -lbmf  -Ivtf 


xxxvxXvxxxxxxx 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


SET  COLOR 

This  is  a  fortran  subroutine  compiled 
separately  to  set  the  color  objects  are 
drawn  in  to  the  bitmap.  The  color  can 
be  set  to  any  of  a  number  of  shade  of 
grey  along  with  black  or  white.  For  this  * 
simulation  the  colors  will  usually 
be  set  to  either  VT  White  or  VT  BLACK. 

This  routine  is  intended  to  be  compiled 
inline  in  a  LISP  environment  to  provide  * 

LISP  graphics  calls.  The  function  is  called  * 

*  as  (setcol  bmd  color)  and  it  will  * 

*  return  a  Mtn  when  complete.  The 
parameters  passed  to  the  fortran  subroutine  * 

m  from  LISP  are:  * 

*  * 

bmd  -  bitmap  descriptor  number 

*  color  =  VT  White 

or 

*  VT  Black 
************************************************** 


* 

Dennis  Poulos  * 

x3403  * 

22  April  1986 

* 


subroutine  setcol(bmd, color) 
integer  bmd. color 
call  bmSetColor(bmd, color) 
return 
end 


* 


compile  as  * 

f77  -c  setcol. f  -Itoolsf  -Ibmf  -Ivtf  * 

* 
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c 


CLEAR  BITMAP 

*  * 

’  This  is  a  fortran  subroutine  compiled 
’  separately  to  clear  a  bitmap  after  it 

*  has  been  drawn  on.  It  will  clear  the  * 

*  region  of  the  bitmap  pointed  to  by  bind  ’ 

*  from  the  position  x  and  y  for  a  width  * 

*  w  and  height  h.  It  will  clear  the  region  * 

*  by  coloring  the  entire  space  in  the  * 

*  color  specified.  These  values  are  passed  * 

*  as  parameters.  This  fortran  subroutine  * 

*  is  intended  to  be  compiled  * 

*  inline  in  a  LISP  environment  to  provide  * 

*  LISP  graphics  calls.  The  function  is  called  * 

*  as  (bmclear  bmd  x  y  w  h  color).  It  will  * 

*  return  a  "t"  when  complete.  The  * 

*  parameters  passed  to  the  fortran  subroutine  * 

*  from  LISP  are:  * 

*  * 

*  bmd  =  bitmap  descriptor  number  * 

*  * 

*  x  =  x  coordinate  of  upper  left  * 

*  corner  of  region  * 

*  * 

*  y  =  y  coordinate  of  upper  left  * 

*  corner  of  region  * 

*  * 

*  w  =  width  of  area  to  be  cleared  * 

*  * 

*  h  =  height  or  area  to  be  cleared  * 

*  * 

*  color  =  color  to  set  the  region  to 

*  to  clear  it.  * 

*  * 

*  all  distances  and  positions  are  measured  * 

*  in  terms  of  pixels  * 

*  * 

************************************************** 

************************************************** 


Dennis  Poulos 

*  x3403 

*  22  April  1986 


subroutine  bmclear) bmd, x,y,w,h, color) 
external  bmClearRegion 
integer  bmd, x,y.w,h. color 
call  bmClearRegion) bmd, x,y,w.h, color) 


SET  THICKNESS 

*  This  is  a  fortran  subroutine  compiled 

*  separately  to  set  the  thickness  of  lines 

*  drawn  onto  the  bitmap.  The  line  thickness  ’ 

*  is  set  in  terms  of  x  number  of  pixels  * 

*  wide.  This  routine  is  intended  to  be  * 

*  compiled  inline  in  a  LISP  environment  to  * 

*  provide  LISP  graphics  calls.  The  function  * 

*  is  called  as  (setthick  bmd  thickness)  * 

*  It  returns  a  "t"  when  complete.  The  * 

*  parameters  passed  to  the  fortran  subroutine  * 

*  from  LISP  are:  * 

*  * 

*  bmd  -  bitmap  descriptor  number  * 

*  * 

*  thickness  =  number  of  pixels  wide  * 

*  the  line  is  desired  to  be  * 

*  * 

************************************************** 

************************************************** 
*  * 

*  Dennis  Poulos  * 

*  x  3403 

*  22  April  1986  * 

*  * 

************************************************** 

subroutine  setthick ( bmd. thickness ) 
integer  bmd. thickness 
call  bmSetThickness( bmd, thickness) 
return 
end 

********  ****************************************** 
*  * 

*  compile  as  * 

*  f77  -c  setthick. f  -ltoolsf  -Ibmf  -lvtf  * 

*  * 

************************************************** 


SET  X  V  POSITION 


*  This  is  a  fortran  subroutine  compiled 
separately  to  position  the  position 

*  pointer  of  a  bitmap.  This  routine  is 

*  intended  to  be  compiled  inline  in  a  LISP  ' 

*  environment  to  provide  LISP  graphics  calls.  * 

*  The  function  is  called  as  * 

*  (bmsetpos  bmd  x  y)  and  it  will  return  a 

’  "t"  when  complete.  The  parameters  passed 

*  to  the  fortran  subroutine  from  LISP  are:  * 

*  * 

*  bmd  =  bitmap  descriptor  number 

x  * 

’  x  =  x  coordinate  of  the  pointer 

x  * 

*  y  =  v  coordinate  of  the  pointer 


*  x 

*  Dennis  Poulos 

*  x  3403 

*  20  April  1986 

X  * 


subroutine  bmsetpos(bmd,x,y) 
integer  bmd.x,y 
call  bmSetPosition(bmd,x,y) 
return 
end 


*  compile  as  ’ 

’  177  -c  bmsetpos. f  -ltoolsf  -lbmf  -lvtf 


1 


S 

A 

»*( 


r********.***** 


SET  LINE  STYLE 


This  is  a  fortran  subroutine  compiled 
separately  to  set  the  style  of  the  line 
to  be  drawl  on  the  bitmap.  This  consists 
of  a  16  bit  mask  as  an  integer  number  and  * 
a  count  that  is  the  number  of  times  each  * 
bit  in  the  mask  is  replacated.  This  routine  * 
is  intended  to  be  compiled  inline  in  a  * 

LISP  environment  to  provide  LISP  graphics  * 
calls.  The  function  is  called  as  * 

(bmsetsty  bmd  mask  count)  and  it  will  * 
return  a  "t"  when  complete.  The  parameters  * 
passed  to  the  fortran  subroutine  from  LISP  * 


*  bmd 


bitmap  descriptor  number 


*  mask 


bit  mask  for  line 


l=  repetition  count  for  bits  in  mask  * 


*  Dennis  Poulos 

*  x  3403 


*  21  April  1986 


subroutine  bmsetsty  (bmd.  mask,count) 
integer  bmd. mask. count 
call  bmSetSty!e(bmd. mask. count) 
return 
end 


*  compile  as  ’ 

*  f77  -c  bmsetsty  . f  -ltoolsf  -Ibmf  -lvtf 
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DRAW  LINE 

*  This  is  a  fortran  subroutine  compiled 

’  separately  to  draw  a  line  from  the  pointer 
'  to  the  x  and  y  position  of  the  end  of  the  * 

*  line.  This  routine  is  intended  to  be 

*  compiled  inline  in  a  LISP  environment  to  * 

*  provide  LISP  graphics  calls.  The  function  * 

*  is  called  as  (bmdrawl  bmd  x  y)  and  it  will  * 

*  return  a  "t"  when  complete.  The  parameters  1 

’  passed  to  the  fortran  subroutine  from  LISP  * 

*  «...  * 
are. 


bmd  =  bitmap  descriptor  number 

* 

x  =  x  coordinate  of  line  end 
v  =  v  coordinate  of  line  end 


*  * 

*  Dennis  Poulos 

*  x  3403 

*  20  April  1986  * 

»  * 


subroutine  bmdrawl(bmd.x,y) 
integer  bmd.x.y 
call  bmPaintLine(bmd.x.y) 
return 
end 

*X***x*****»********«********S**X**X*************X 
*  * 

*  compile  as  * 

*  f77  -c  bmdrawl. f  -Itoolsf  -Ibmf  -lvtf 

x  x 

****xxx*m»*«*x***»xxx***«****x**xX*x*xx*******x**x 
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DISPLAY  BIT  MAP 


This  is  a  fortran  subroutine  compiled 
separately  to  display  a  bitmap  in  a  window. 
This  routine  is  intended  to  be  compiled 
inline  in  a  LISP  environment  to  provide 
LISP  graphics  calls.  The  function  is  called  ’ 
"bmdisp".  From  LISP  it  is  called  as  * 

(bmdisp  fd  bmd  mode  sx  sy  dx  dy  w  h  color) 
and  it  will  return  0  which  equals  the  * 
copyraster  function.  The  parameters  passed 
to  the  fortran  function  from  LISP  are:  * 


fd  =  file  descriptor  number 


bmd  =  bitmap  descriptor  number 


type  --  type  of  bmdisplay  (copyraster 
toggleraster.  paintraster)  * 


x  coordinate  of  upper  left 
corner  of  portion  of  bitmap 


y  coordinate  of  upper  left 
corner  of  portion  of  bitmap 


x  coordinate  of  upper  left 
corner  of  destination  window 


=  y  coordinate  of  upper  left 
corner  of  destination  window 


width  of  region  (source  and 
display)  * 


height  of  region  (source  and 
display)  * 


color  (used  if  using  paintraster,  * 
otherwise  set  to  0) 


Dennis  Poulos 
x  3403 

19  April  19t>6 


integer  function  bmdisp(fd. type. bmd.sx.sy  .dx.dy  .w.h. color) 
external  bmDisplayBitmap 

integer  bmDisplavBitmap.fd.  type,  bmd.sx.sy. dx.dy.w.h.  color 
bmdisp  =  bmDisplayBumap(fd.type.bmd.sx.sy.dx.dy  .w  ,h. color) 

return 

end 

****X*****X*X****«S*XX***X**«X*X3|C*X********«**XX** 


c  *  compile  as  * 

c  *  f77  -c  bmdisp. f  -ltooisf  -lbmf  -lvtf  * 


c 


K  *  *  *  *  *  * 


i*r**K«* 
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CLOSE- WINDOW 


*  This  is  a  fortran  subroutine  compiled 

*  separately  to  close  a  window  on  the  ISI 

*  Graphics  Machines.  It  is  intended  to  be 

’  compiled  inline  in  a  LISP  environment  to 

*  provide  LISP  graphics  calls.  The  function 

*  is  called  "clw".  From  LISP  it  is  called 

*  as  (clw  fd)  and  it  will  return  a  "t"  when 

*  complete.  The  parameters  passed  to  the 

*  fortran  function  from  LISP  are:  * 


fd  =  file  descriptor  of  the  desired  * 
window  * 


Dennis  Poulos 
x  3403 

15  April  1986 


subroutine  clw(fd) 
integer  fd 

call  CloseWindow(fd) 


compile  as 

f77  -c  closew.f  -ltools  -lbmf  -lvtf 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


xxxxxXXMXXXXXxxxxmxxxxx*****************-********** 


BITMAP  DEALLOCATE 

*  This  is  a  fortran  subroutine  compiled 

*  separately  to  deallocate  space  in  memory 

*  for  a  bitmap.  A  bit  map  is  maintained  in  * 

*  memory  and  can  be  displayed  in  a  window  * 

*  as  desired.  This  routine  is  intended  to  be  * 

*  compiled  inline  in  a  LISP  environment  to  * 

*  provide  LISP  graphics  calls.  The  function  * 

*  is  called  ’'bmdeal".  From  LISP  it  is  called  * 

*  as  (bmdeal  bind)  and  it  will  return  a  "t"  * 

*  when  complete.  The  parameters  passed  to  * 

*  the  fortran  subroutine  from  LISP  are:  * 

*  X 

*  bmd  =  bitmap  descriptor  * 

*  X 

****x**********«************x*******x**x«********* 


X**X**X***********«*X***«***************X********* 
X  X 

’  Dennis  Poulos  * 

*  x  3403  * 

*  19  April  1986 

X  X 


subroutine  bmdeal(bmd) 
integer  bmd 

call  bmDeallocate(bmd) 
return 
end 


compile  as  * 

177  -c  bmdeal. f  -ltoolsf  -lbmf  -lvtf  * 

* 
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APPENDIX  C  -  LSCAN 


* 

This  is  a  lisp  program  to  load  and  execute 
graphics  on  the  IS1  machines  and  to  do  image  * 

processing.  * 

* 

******************************************************* 


******************************************************* 

* 

GRAPHICS  ROUTINES 

* 

The  following  "cfasl"  lisp  functions  compile  * 
inline  the  ISI  graphics  routines.  They  use  * 
fortran  libraries  in  the  ISI  machines.  The  * 
first  term  after  the  "cfasl"  is  the  file  name  * 
of  the  object  code  produced  by  the  fortran  * 
compiler.  The  second  term  is  the  entry  point  * 
into  the  subroutine.  The  third  term  is  the  name  * 
given  to  the  LISP  function  created  by  the  * 

"cfasl".  The  forth  term  is  type  of  routine  it  * 
is  and  the  last  term  are  the  libraries  the  * 

compiler  needs  to  link  to.  * 

* 

***»**••***•******•***•»*****••****«*•****•**••*****«** 


(cfasl  ’tofile.o  'tofile  ’tofile  nil  "-Ivt  -lm  -Ibm") 

(cfasl  ’openw.o  ’_opw_  ’opw  "integer-function"  "-ltoolsf  -lbmf  -lvtf") 

(cfasl  'setldisc.o  ’_slined_  ’slined  "integer-function"  "-ltoolsf  -lbmf  -lvtP') 
(cfasl  ’bitmal.o  ’_bitmapa_  ’bitmaps  "integer-function"  "-ltoolsf  -lbmf  -lvtf') 
(cfasl  ’bmadres.o  ’_bmadres_  ’bmadres  "subroutine"  "-ltoolsf -lbmf  -lvtf") 
(cfasl  "writfile.o  ’_writfile_  ’writ file  "integer-function"  "-ltoolsf  -lbmf -lvtP') 
(cfasi  'setbacc.o  ’_setbacc_  ’setbacc  "subroutine"  "-ltoolsf  -lbmf  -lvtf') 

(cfasl  ’setcol.o  '_setcol_  ’setcol  "subroutine"  "-ltoolsf  -lbmf  -lvtP’) 

(cfasl  ’bmclear.o  ’_bmclear_  ’bmclear  "subroutine"  "-ltoolsf  -lbmf  -lvtf') 
(cfasl  ’setthick.o  ’_setthick_  ’setthick  "subroutine"  "-ltoolsf  -lbmf  -lvtf') 

(cfasl  ’bmsetpos.o  ’_bmsetpos_  ’bmsetpos  "subroutine"  "-ltoolsf  -lbmf  -lvtP') 
(cfasl  ’bmsetsty.o  ’_bmsetsty_  ’bmsetsty  "subroutine"  "-ltoolsf  -lbmf  -lvtf") 
(cfasl  ’bmdrawl.o  ’_bmdrawl_  ’bmdrawl’  "subroutine"  "-ltoolsf  -lbmf  -lvtf") 
(cfasl  ’bmdisp.o  ’_bmdisp_  ’bmdisp  "integer-function"  "-ltoolsf  -lbmf  -lvtP' ) 
(cfasl  'closew.o  ’_clw_  ’clw  "subroutine"  "-ltoolsf  -lbmf  -IvtF') 

(cfasl  ’bmdeal.o  ’_bmdeal_  ’bmdeal  "subroutine"  "-ltoolsf  -lbmf  -lvtP’) 

******************************************************** 

*  * 

CONSTANTS 
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The  following  are  constants  used  throughout 
the  simulation. 


(setq  wwidth  576) 

; Window  Wridth 

(setq  wheight  448) 

•  Window  Height 

(setq  filenum  40) 

:a  file  counter  to  create 
;new  files  to  copy  the  bit 
•.maps  to 

(setq  filebase  filenum) 

;base  number  for  the  files 
;set  to  file  number 

(setq  numruns  5) 

;the  number  of  desired  runs 

(setq  wlx  35) 

;X  coordinate  of  the  upper 
deft  corner  of  window  1 

(setq  wly  50) 

;Y  coordinate  of  the  upper 
deft  corner  of  window  1 

(setq  xsep  54) 

dateral  separation  between 
;windows 

(setq  ysep  52) 

ivertical  separation  between 
:windows 

(setq  pie  3.14159) 

;P' 

(setq  slope  limit  2.0) 

dimit  between  flat  terrain 
;and  a  sloping  terrain 

(setq  curv  threshold  0.05) 

threshold  for  measurable 
;  curvature 

(setq  unsafe  slope  limit  30.0)  threshold  between  safe  and 

;  unsafe  slope  in  degrees 

(setq  xOinit  287) 

;terrain  bitmap  xO  coordinate 

(setq  yOinit  133) 

;terrain  bitmap  yO  coordinate 

(setq  xOcminit  77) 

classification  map  xO 
coordinate 

(setq  yOcminit  13) 

classification  map  yO 
coordinate 

(setq  peak  char  ’P) 

;peak  symbol  for 
classification  map 

(setq  pit  char  ’L ) 

;pit  symbol  for 
classification  map 

(setq  ridge  char  ’  “) 

;ridge  symbol  for 
classification  map 

(setq  ravine  char  v) 

;ravine  symbol  for 
classification  map 

(setq  saddle  char  'S) 

isaddle  symbol  for 
classification  map 

(setq  pass  char  ’Y) 

;pass  symbol  for  the 
classification  map 

(setq  flat  char  ’-) 

;flat  terrain  symbol  for 
classification  map 

(setq  safe  slope  char  ’/) 

.safe  slope  symbol  for 
classification  map 

(setq  unsafe  slope  char  ’U) 

;unsafe  slope  symbol  for 

€ 


(setq  edgechar  'E) 

(setq  numacross  29) 

(setq  numdown  29) 
(setq  bmdx  9) 

(setq  bmdy  3) 

(setq  windows  nil) 

(setq  rclassmap  nil) 

(setq  nclassmap  nil) 

(array  real_terrain  t  (1-r- 

(array  noisy  terrain  t  (1 

(setq  noiselimit  1) 


(setq  full  line  65535) 
(setq  VT  ABSOLUTE  1) 
(setq  VT  RELATIVE  0) 
(setq  COPYRASTER  0) 
(setq  TWSDISC  5) 

(setq  VT  BLACK  0) 
(setq  VT  WHITE  I) 


:classification  map 

■.edge  symbol  for 
:classification  map 

:0->numacross  terrain 
:cells  across 

:0->numdown  terrain 
:cells  down 

;bitmap  delta  x  per 
:terrain  cell 

;  bitmap  delta  y  per 
;terrain  cell 

;sets  list  of  windows 
:equal  to  the  nil  set 
;sets  the  list  of 
:real  _terrain  classification 
:  map  pixels  to  nil 

;sets  the  list  of 
: noisy  terrain  classification 
.map  pixels  to  nil 
numacross) (1-^  numdown)  3) 

:30  x  30  x  3  array  of  map 
:x,y.z  coordinates 
-  numacross)(l—  numdown)  3) 

:30  x  30  x  3  array  of 
:noisy  x,y,z  map 
coordinates 

;e)evation  data  (z)  is 
.subject  to  ■+/-  the 
:noise  limit  in 
;noise  measured  in  pixels 
;mask  for  a  full  line 

;addressing  mode 
;relative  addressing  mode 
.-bitmap  display  mode 
dine  discipline  mode 
; black  color 
;white  color 


USER  DEFINED  LISP  FUNCTIONS 

x  * 

*  The  following  are  the  user  defined  functions  * 

*  that  do  the  actual  work  in  this  simulation.  * 

*  * 

******************************************************** 

**  "copy  _bm  tofile"  copys  the  bitmaps  of  the  terrain 
**  and  the  symbolic  maps  to  files  so  that  they  may  be 
**  printed  of  the  laser  printer 

(defun  copy_bm_tofile  (bmd) 
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(setq  filenum  (1-  filenum)) 

(tofile  bmd  (getpname  (concat  'figure,  filenum)))) 

”  "def  window"  create?  a  list  of  the  open  windows 
**  and  attaches  attributes  to  each  window 

(defun  def  window  (winnum  xcoord  ycoord  width  height  name) 
(putprop  winnum  xcoord  'xcoord) 

(putprop  winnum  ycoord  'ycoord) 

(putprop  winnum  width  ’width) 

(putprop  winnum  height  ’height) 

(putprop  winnum  name  ’name) 

(setq  windows  (appendl  windows  winnum))) 

**  "xtrfunc"  defines  the  transfer  function  for 
**  the  terrain  map  x  y  z  coordinates  to  be 
**  mapped  into  the  window  x  y  coordinates 

(defun  xtrfunc  (x  y) 

(—  xOinit  (-  (*  bmdx  x)(*  bmdx  y))}) 

*’  "ytrfunc"  defines  the  transfer  function  for 
**  the  terrain  map  x  y  z  coordinates  to  be 
**  mapped  into  the  window  x  y  coordinates 

(defun  ytrfunc  (x  y  z) 

(-  (-  yOinit  (*  bmdy  x)(*  bmdy  y))  z)) 

**  "arrayx"  recovers  the  x  value  of  a  terrain  array 
**  for  a  given  pixel  pi  pj 

(defun  arrayx  (pi  pj  array  name) 

(array  call  t  array  name  pi  pj  0)) 

**  "arrayy"  recovers  the  y  value  of  a  terrain  array 
**  for  a  given  pixel  pi  pj 

(defun  arrayy  (pi  pj  arrayname) 

(arraycall  t  arrayname  pi  pj  1)) 

**  "arrayz"  recovers  the  z  value  of  a  terrain  array 
*'  for  a  given  pixel  pi  pj 

(defun  arrayz  (pi  pj  arrayname) 

(arraycall  t  arrayname  pi  pj  2)) 

**  "paint  saddle"  paints  the  character  "S"  on  the 
**  symbolic  map  to  represent  a  saddle. 

(defun  paint  saddle  (bmd  cmx  cmy) 

(bmsetpos  bmd  (*  10  cmx)  (+  2  cmy)) 

(bmdrawl  bmd  -7  0) 
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(bmdrawl  bind  0  5) 

(bmdrawl  bind  7  0) 

(bmdrawl  bmd  0  4) 

(bmdrawl  bmd  -7  0)) 

**  "paint_pass"  paints  the  character  "P"  on  the 
*’  symbolic  map  to  represent  a  pass. 

(defun  paint _pass  (bmd  cmx  cmy) 

(bmsetpos  bmd  (-<-  3  cmx)  (i-  11  cmy)) 

(bmdrawl  bmd  0  -9) 

(bmdrawl  bmd  5  0) 

(bmdrawl  bmd  2  2) 

(bmdrawl  bmd  0  1) 

(bmdrawl  bmd  -2  2) 

(bmdrawl  bmd  -4  0)) 

;  **  "paint  safe  slope"  paints  the  character  /  on  the 

**  symbolic  map  to  represent  a  safe  slope. 

(defun  pain  t_safe  slope  (bmd  cmx  cmy) 

(bmsetpos  bmd  (-  cmx  2)  (+  cmy  11)) 

(bmdrawl  bmd  9  -9)) 

**  "paint  edge"  paints  the  character  "E"  on  the 
**  symbolic  map  to  represent  an  edge. 

(defun  paint  edge  (bmd  cmx  cmy) 

(bmsetpos  bmd  ( —  cmx  10)  (+  cmy  2)) 

(bmdrawl  bmd  -7  0) 

(bmdrawl  bmd  0  5) 

(bmdrawl  bmd  4  0) 

(bmdrawl  bmd  -4  0) 

(bmdrawl  bmd  0  4) 

(bmdrawl  bmd  7  0)) 

**  "paint  unsafe  slope"  paints  the  character  "U"  on  the 
**  symbolic  map  to  represent  an  unsafeslope. 

(defun  paint  unsafe  slope  (bmd  cmx  cmy) 

(bmsetpos  bmd  (f  cmx  3)  (+  cmy  2)) 

(bmdrawl  bmd  0  6) 

(bmdrawl  bmd  1  2) 

(bmdrawl  bmd  2  1) 

(bmdrawl  bmd  1  0) 

(bmdrawl  bmd  2-1) 

(bmdrawl  bmd  1  -2) 

(bmdrawl  bmd  0  -6)) 

**  "paint  peak"  paints  the  character  for  a  peak  symbol 
**  on  the  symbolic  map  to  represent  a  peak. 
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(defun  paint  _peak  (bmd  cmx  cmy) 

(bmsetpos  bmd  (—  cmx  2)  (—  cmy  6)) 

(bmdrawl  bmd  4  -4 1 
(bmdrawl  bmd  0  9) 

(bmdrawl  bmd  0  -9) 

(bmdrawl  bmd  4  4)) 

;  **  "paintpit"  paints  the  character  for  a  pit  symbol 

**  on  the  symbolic  map  to  represent  a  pit. 

(defun  paint  pit  (bmd  cmx  cmy) 

(bmsetpos  bmd  (—  cmx  2)  (-t-  cmy  7)) 

(bmdrawl  bmd  4  4) 

(bmdrawl  bmd  0  -9) 

(bmdrawl  bmd  0  9) 

(bmdrawl  bmd  4-4)) 

**  "paint_ravine"  paints  the  character  "V"  on  the  symbolic 
**  map  to  represent  a  ravine  on  the  symbolic  map. 

(defun  paintravine  (bmd  cmx  cmy) 

(bmsetpos  bmd  (—  cmx  2)  (•+  cmy  2)) 

(bmdrawl  bmd  4  9) 

(bmdrawl  bmd  1  0) 

(bmdrawl  bmd  4  -9)) 

;  **  "paint  ridge"  paints  the  symbol  on  the  symbolic  map 

**  to  represent  a  ridge  on  the  symbolic  map. 

t 

(defun  paint  ridge  (bmd  cmx  cmy) 

(bmsetpos  bmd  (—  cmx  2)  (-r  cmy  11)) 

(bmdrawl  bmd  4  -9) 

(bmdrawl  bmd  1  0) 

(bmdrawl  bmd  4  9)) 

**  "paint  level"  paints  the  character  on  the  symbolic 
**  map  to  represent  a  level  terrain  on  the  symbolic  map. 

(defun  paint  flat  (bmd  cmx  cmy) 

(bmsetpos  bmd  (—  cmx  3)  (  —  cmy  7)) 

(bmdrawl  bmd  7  0)) 

;  **  "make  real  terrain"  is  a  function  that  fills  the 

;  **  array  real_terrain  with  x,y,z  coordinates  for  the 

**  synthetic  terrain  used  in  this  model.  At  this  time 
**  it  does  not  create  any  elevation  (z).  This  will  be 

;  **  added  later. 

(defun  make  real  terrain  () 

(prog  (x  y  z  pi  pj  rl  r2  r3  r4  r5  counter) 

(setq  x  0) 

(setq  y  0) 
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(setq  pj  0) 

(setq  counter  OJ 
(setq  r5  ( —  2  (random  5))) 
loop 

(store  (real_terrain  pi  pj  0)  x) 

(store  (realterrain  pi  pj  1)  y) 

(store  (real  terrain  pi  pj  2)  z ) 

(cond  ((  and  (=  pi  numacross  )(=  pj  namdown  ) )(go  loopl)) 
((=  pi  numacross  ) 

(setq  pi  0) 

(setq  x  0) 

(setq  pj  (1+  pj)) 

(setq  y  (1-  y)) 

(go  loop) ) 

(t  (setq  pi  (1- pi)) 

(setq  x  (1+  x)) 

(go  loop))) 
loopl 

(cond  ((=  counter  r5)(return))) 

(setq  pi  0) 

(setq  pj  0) 

(setq  rl  (random  31)) 

(setq  r2  (random  31)) 

(setq  r3  (quotient  (float  (-*-  5  (random  11)))  10.0  )) 

(setq  r4  (float  (random  51))) 
loop2 

(setq  z  ( array z  pi  pj  ’real  terrain)) 

(cond  ((evenp  counter) 

(setq  z  (t-  z  (fix  (times  r4  (  cos 

(times  r3  2.0  pie(quotient  (sqrt 
(add  (float  (expt  (-  rl  pi)  2)) 

(float  (expt  (-  r2  pj)  2))))  30  )))))))) 

(t 

(setq  z  (-i-  z  (fix  (times  r4  (  sin 

(times  rS  2.0  pie 
(quotient  (sqrt  (add  (float 

(expt  (-  rl  pi)  2)) 

(float  (expt  (-  r2  pj)  2))))  30  ))))))))) 

(store  (real  terrain  pi  pj  2)  z) 

(cond  ((and  (=  pi  numacross)(=  pj  namdown)) 

(setq  counter  (l-<-  counter)) 

(go  loopl)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1*  pj)) 

(go  loop2)) 

(t 

(setq  pi  (1-  pi) ) 

(go  loop2))))) 


*  * 


"makenoisy  terrain"  takes  the  arrav  real  terrain 
**  and  adds  a  random  noise  to  the  elevation  data.  The 
**  noise  is  —  the  constant  noise  limit 

(defun  make  noisy  terrain  () 

(prog  (pi  pj) 

(setq  pi  0) 

(setq  pj  0) 
loop 

(store  (noisy  terrain  pi  pj  0)(&rrayx  pi  pj  ’realterrain)) 
(store  (noisy  terrain  pi  pj  l)(arrayy  pi  pj  ’real  terrain)) 
(store  (noisy  terrain  pi  pj  2) 

(—  (arrayz  pi  pj  ’real_terrain) 

(-  noise  limit  (random  (1-r-  (*  2  noise_limit )))))) 
(cond  ((  and  (=  pi  numacross)(=  pj  numdown)) 

(return)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1+  pj)) 

(go  loop)) 

(t 

(setq  pi  (1*  pi)) 

(go  loop))))) 

**  "initialize"  opens  the  required  four  windows 
**  and  attaches  the  fd  (file  descriptor)  number 
**  and  the  bmd  (bitmap  descriptor)  number  to  the 
**  window-list  as  attributes  of  each  window.  It  also 
**  sets  the  line  discipline  of  each  window  to  the 
**  graphics  mode  vice  the  standard  TTY  mode. 

(defun  initialize  () 

(prog  (windowlist  x) 

(setq  windowlist  windows) 
loop 

(setq  x  (car  windowlist)) 

(setq  windowlist  (cdr  windowlist)) 

(putprop  x  (opw  (get  x  ’xcoord) 

(get  x  'y coord) 

(get  x  ’width) 

(get  x  ’height) 

(get  x  ’name)) 

’fdnum) 

(putprop  x  (bitmaps  (get  x  ’width)(get  x  ’height))  ’bmd) 
(bmadres  (get  x  ’bmd)  VT_ABSOLUTE) 

(setbacc  (get  x  ’bmd)  VT  WHITE) 

(setcol  (get  x  ’bmd)  VT  BLACK) 

(bmsetsty  (get  x  ’bmd)  full _line  2) 

(bmciear  (get  x  ’bmd)  0  0  wwidth  wheight  VT_WHITE) 
(slined  (get  x  ’fdnum)  TWSDISC) 

(cond  ((  not(  null  windowlist))  (go  loop))))) 


*  * 


"draw  terrain  base"  draws  the  base  on  which  the 
synthetic  terrain  is  overlaid 

(defun  draw  terrain  base  (bmd) 

(prog  (xl  xm  yu  xr  ym  yl  xr) 

(setq  xm  xOinit) 

(setq  yu  yOinit) 

(setq  xr  (—  xOinit  (’  numacross  bmdx))) 

(setq  ym  (—  yOinit  (*  numdown  bmdy  ))) 

(setq  yl  (—  yOinit  (*  2  numdown  bmdy))) 

(setq  xl  (-  xOinit  (*  numacross  bmdx))) 

(bmsetpos  bmd  xm  yu) 

(setthick  bmd  2) 

(bmdrawl  bmd  xr  ym) 

(bmdrawl  bmd  xm  yl) 

(bmdrawl  bmd  xl  ym) 

(bmdrawl  bmd  xm  yu) 

(setthick  bmd  1))) 

’*  "draw  terrain"  takes  the  x  y  z  coordinates  of 
**  the  terrain  held  in  the  terrain  arrays  and 
**  plots  them 

(defun  draw  terrain  (arrayname  bmd  fd) 

(prog  (pi  pj  x  y  z  bmx  bmy) 

(setq  pi  0) 

(setq  pj  0) 

(draw  terrain  base  bmd) 
loop 

(setq  x  (arrayx  pi  pj  arrayname)) 

(setq  y  (arrayy  pi  pj  arrayname)) 

(setq  z  (arrayz  pi  pj  arrayname)) 

(setq  bmx  (xtrfunc  x  y)) 

(setq  bmy  (ytrfunc  x  y  z)) 

(bmsetpos  bmd  bmx  bmy) 
loopl 

(cond((and(=  pi  numacross)(=  pj  numdown)} 

(setthick  bmd  2) 

(bmdrawl  bmd  bmx  (+  bmy  z) ) 

(bmsetpos  bmd  bmx  bmy) 

(setthick  bmd  1) 

(bmdisp  fd  COPYRASTER  bmd  0  0  0  0  wwidth  wheight  0) 
(return)) 

((=  pj  numdown) 

(setthick  bmd  2) 

(bmdrawl  bmd  bmx  (-*-  bmy  z)) 

(bmsetpos  bmd  bmx  bmy) 

(setthick  bmd  1) 

(setq  pi  (1*  pi)) 

(setq  x  (arrayx  pi  pj  arrayname)) 

(setq  y  (arrayy  pi  pj  arrayname)) 

(setq  z  (arrayz  pi  pj  arrayname)) 
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(seiq  bmy  (vtrfunc  xy  z)) 

(bmdrawl  bmd  bmx  bmy) 

(go  loopl)) 

((=  pi  numacrossj 
(setthick  bmd  2) 

(bmdrawl  bmd  bmx  (—  bmy  z)) 

(bmsetpos  bmd  bmx  bmy) 

(setthick  bmd  1) 

(seiq  pj  (1-  pj)) 

(setq  x  (arrayx  pi  pj  array  name)) 

(setq  y  (arrayy  pi  pj  arrayname)) 

(setq  z  (arrayz  pi  pj  arrayname)) 

(bmdrawl  bmd  (xtrfunc  x  y)(ytrfunc  x  y  z)) 

(setq  pi  0) 

(go  loop)) 

((or  (=  pj  0)(=  pi  0)) 

(setthick  bmd  2) 

(bmdrawl  bmd  bmx  (—  bmy  z)) 

(bmsetpos  bmd  bmx  bmy) 

(setthick  bmd  1) 

(setq  pi  ( 1  —  pi) ) 

(setq  x  (arrayx  (1-  pi) ( 1  -r-  pj)  arrayname)) 

(setq  y  (arrayy  (1-  pi ) ( 1  pj)  arrayname)) 

(setq  z  (arrayz  (1-  pi) ( 1  —  pj)  arrayname)) 

(bmdrawl  bmd  (xtrfunc  x  y)(ytrfunc  x  y  z)) 
(bmsetpos  bmd  bmx  bmy) 

(setq  x  (arrayx  pi  pj  arrayname)) 

(setq  y  (arrayy  pi  pj  arrayname)) 

(setq  z  (arrayz  pi  pj  arrayname)) 

(setq  bmx  (xtrfunc  x  y ) ) 

(setq  bmy  (ytrfunc  x  y  z)) 

(bmdrawl  bmd  bmx  bmy) 

(go  loopl)) 

(t 

(setq  pi  (1-  pi)) 

(setq  x  (arrayx  (1-  pi)(l-*-  pj)  arrayname)) 

(setq  y  (arrayy  (1-  pi)(l~  pj)  arrayname)) 

(setq  z  (arrayz  (1-  pi ) ( 1  — -  pj)  .arrayname)) 

(bmdrawl  bmd  (xtrfunc  x  y) (ytrfunc  x  y  z)) 
(bmsetpos  bmd  bmx  bmy) 

(setq  x  (arrayx  pi  pj  arrayname)) 

(setq  y  (arrayy  pi  pj  arrayname)) 

(setq  z  (arrayz  pi  pj  arrayname)) 

(setq  bmx  (xtrfunc  x  y ) ) 

(setq  bmy  (ytrfunc  xy  z)) 

(bmdrawl  bmd  bmx  bmy) 

(go  loopl))))) 

**  "drawclass  map  base"  draws  the  base  in  which  the 

**  classification  map  is  displayed 
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no-0171  851 

UNCLASSIFIED 


RANGE  IMAGE  PROCESSING  FOR  LOCAL  NAVIGATION  OF  AN  1/1 

AUTONOMOUS  LAND  VEHICLES)  NAVAL  POSTGRADUATE  SCHOOL 
MONTEREV  CA  D  D  POULOS  SEP  86 

F/G  6/4  NL 


(defun  draw_class_map_base  (bmd) 

(prog  () 

(setthick  bmd  2) 

(bmsetpos  bmd  (•  xOcminit  5)  (-  yOcminic  5)) 

(bmadres  bmd  VT  RELATIVE)' 

(bmdrawl  bmd  430  0) 

(bmdrawl  bmd  0  430) 

(bmdrawl  bmd  -430  0) 

(bmdrawl  bmd  0  -430) 

(setthick  bmd  1))) 

**  "makerclassmap"  builds  the  list  of 
**  pixels  associated  with  a  classification 
**  map.  It  will  determine  if  a  pixel  needs  to  be 
*  '  classified  or  if  it  is  a  edge  pixel. 

**  If  it  needs  to  be  classified  it  will 
**  invoke  the  function  "classify". 

(defun  make  rdass  map  () 

(prog  (pi  pj  x  y  z  ) 

(setq  pi  0) 

(setq  pj  0) 

(setq  arrayname  ’real_terrain) 
loop 

(setq  pixname  (concat  ’pixel-  (+  pi  (*  30  pj)))) 

(putprop  pixname  (arrayx  pi  pj  arrayname)  ’x) 

(putprop  pixname  (arrayy  pi  pj  arrayname)  ’y) 

(putprop  pixname  (quotient  (float  (arrays  pi  pj  arrayname) 

)  10.0)  ’*) 

(cond  ((or  (=  pi  0)(=  pi  numacross)(=  pj  0)(=  pj  numdown)) 
(putprop  pixname  ’edge  ’pri  class) 

(putprop  pixname  edgechar  ’letter)) 

(t 

(classify  pixname  arrayname  pi  pj))) 

(setf  rclass  map  (appendl  rclass  map  pixname)) 

(cond  ((and  (=  pi  numacross)(=  pj  numdown)) 

(return)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1+  pj)) 

(go  loop)) 

(t 

(setq  pi  (1—  pi)) 

(go  loop))))) 

**  "makenclassmap"  builds  the  list  of  pixels  associated 
**  with  a  classification  map.  It  will  determine  if  a 
**  pixel  needs  to  be  classified  or  if  it  is  a  edge  pixel. 

**  If  it  needs  to  be  classified  it  will  invoke  the 
**  function  "classify". 

(defun  make  nclass  map  () 
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(prog  (pi  pj  x  y  z  ) 

(setq  pi  0) 

(setq  pj  0) 

(setq  arrayname  'noisy  _terrain) 
loop 

(setq  pixname  (concat  ’pixel*  (—  pi  (*  30  pj)))) 

(putprop  pixnaine  (arrayx  pi  pj  arrayname)  ’x) 

(putprop  pixname  (arrayy  pi  pj  arrayname)  ’y) 

(putprop  pixname  (quotient  (float  (arrayz  pi  pj  arrayname) 

)  10.0)  ’z) 

(cond  ((or  (=  pi  0)(=  pi  numacross)(=  pj  0)(=  pj  numdown)) 
(putprop  pixname  ’edge  ’pri  class) 

(putprop  pixname  edgechar  ’letter)) 

(t 

(classify  pixname  arrayname  pi  pj))) 

(setf  nc lass  map  (appendl  nclassmap  pixname)) 

(cond  ((and  (  =  pi  numacross)(=  pj  numdown)) 

(return)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1+  pj)) 

(go  loop)) 

(t 

(setq  pi  (1+  pi)) 

(go  loop))))) 

;  **  "classify"  does  the  real  work  of  generating  the 

**  max  and  min  curvature  and  the  slope  of  each 
;  **  pixel  and  deciding  what  it  means. 

(defun  classify  (pixname  arrayname  pi  pj) 

(prog  (i  j  x  y  z  k2  kS  k4  k5  k6  elev  b  c  El  E2  temp  slope) 

(setq  i  (1*  pi)) 

(setq  j  (1-  pj)) 

(setq  x  0) 

(setq  y  0) 

(setq  z  0) 

(setq  k2  0) 

(setq  k3  0) 

(setq  k4  0) 

(setq  k5  0) 

(setq  k6  0) 

(setq  elev  0) 
loopl 

(setq  x  (diff  (float  (arrayx  i  j  arrayname)) 

(float  (arrayx  pi  pj  arrayname)))) 

(setq  y  (diff  (float  (arrayy  i  j  arrayname)) 

(float  (arrayy  pi  pj  arrayname)))) 

(setq  z  (quotient 

(diff  (float  (arrayz  i  j  arrayname)) 

(float  (arrayz  pi  pj  arrayname))) 

10.0)) 
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(setq  k2  (plus  k2  (times  x  z) )) 

(setq  k3  (plus  k3  (times  y  z ) ) ) 

(setq  k4  (plus  k4  (times  x  x  z) ) ) 

(setq  ko  (plus  k5  (times  x  y  z ) ) ) 

(setq  k6  (plus  k6  (times  y  y  z) ) ) 

(setq  elev  (plus  elev  z)) 

(cond  ((=  i  (1+  pi)) 

(setq  i  (1-  pi)) 

(setq  j  (1+  j)) 

(go  loopl)) 

((and  (nequal  i  (1+  pi))  (nequal  j  (1+  pj) ) ) 

(setq  i  (1+  i)) 

(go  loopl))) 

(setq  k2  (quotient  k2  6.0)) 

(setq  k3  (quotient  kS  6.0)) 

(setq  k4  (diff  (quotient  k4  2.0) 

(quotient  elev  3.0))) 

(setq  k5  (quotient  k5  4.0)) 

(setq  k6  (diff  (quotient  k6  2.0) 

(quotient  elev  3.0))) 

(setq  b  (minus  (sum  (times  2.0  k4)  (times  2.0  k6)))) 

(setq  c  (diff  (times  4.0  k4  k6)(times  kS  k5))) 

(setq  El  (quotient  (sum  (minus  b) 

(sqrt  (diff  (times  b  b) (times  4.0  c)))) 

2.0)) 

(setq  E2  (quotient  (diff  (minus  b) 

(sqrt  (diff  (times  b  b)(times  4.0  e)))) 

2.0)) 

(cond  ((greaterp  (abs  E2)(abs  El)) 

(setq  temp  El) 

(setq  El  E2) 

(setq  E2  temp))) 

(setq  slope  (quotient  (times  360.0 

(atan  (sqrt  (sum  (times  k2  k2) 

(times  k3  k3)))  1.0)) 

2.0  pie)) 

(putprop  p ix name  slope  ’slope) 

(putprop  pixname  El  ’maxcurvature) 

(putprop  pixname  E2  ’min  curvature) 

FIRST  RULE 

-  IF  SLOPE  IS  APPROX  0  -  THE  PRIMARY  CLASSIFICATION 
OF  THE  PIXEL  IS  LEVEL  TERRAIN 
(cond  ((<  slope  sicpe  limit) 

(putprop  pixname  ’level  ’pri  class) 

-AND  IF  MAX  AND  MIN  CURVATURE  IS  NEGATIVE 
-PEAK  IS  SECONDARY  CLASSIFICATION 
-PEAK  IS  LETTER 

(cond  ((and  (<  El  (minus  curv  threshold)) 

(<  E2  (minus  curv_threshold))) 

(putprop  pixname  ’peak  ’sec  class) 

(putprop  pixname  peak  char  ’letter)) 

;OR  MAX  AND  MIN  CURVATURE  IS  POSITIVE 


:-PIT  IS  SECONDARY  CLASSIFICATION 
:-PIT  IS  'LETTER 

((and  (>  El  curv_threshold) 

(>  E2  curv_threshold)) 

(putprop  pixname  'pit  ’sec  class) 

(putprop  pixname  pit  char  ’letter)) 

;OR  MAX  CURVATURE  IS  NEGATIVE  WHILE  MIN  CURVATURE  IS  APPROX  0 
RIDGE  IS  SECONDARY  CLASSIFICATION 
:-RIDGE  IS  'LETTER 

((and  (<  El  (minus  curvthreshold)) 

(<  (abs  E2)  curv_threshold)) 

(putprop  pixname  'ridge  ’sec_elass) 

(putprop  pixname  ridge  char  'letter)) 

:OR  MAX  CURVATURE  IS  POSITIVE  WHILE  MIN  CURVATURE  IS  APPROX  0 
:-RAVINE  IS  SECONDARY  CLASSIFICATION 
>RAVINE  IS  LETTER 

((and  (>  El  curv  threshold) 

(<  (abs  E2)  curv  threshold)) 

(putprop  pixname  'ravine  ’sec  class) 

(putprop  pixname  ravine  char  ’letter)) 

OR  IF  THE  MAX  CURVATURE  IS  NEGATIVE  WHILE  THE  MIN  CURVATURE 
:  POSITIVE 

>SADDLE  IS  SECONDARY  CLASSIFICATION 
:-SADDLE  IS  'LETTER 

((and  (<  El  (minus  curv  threshold)) 

(>  E2  curv  threshold)) 

(putprop  pixname  'saddle  'sec  class) 

(putprop  pixname  saddle  char  ’letter)) 

OR  IF  THE  MAX  CURVATURE  IS  POSITIVE  WHILE  THE  MIN  CURVATURE 
:IS  NEGATIVE 

:-PASS  IS  SECONDARY  CLASSIFICATION 
-PASS  IS  'LETTER 

((and  (>  El  curv  threshold) 

(<  E2  (minus  curv  threshold))) 

(putprop  pixname  'pass  ’sec_class) 

(putprop  pixname  pass  char  'letter)) 

:EVERYTHING  ELSE  IS 

-FLAT  TERRAIN  IS  SECONDARY  CLASSIFICATION 
;-FLAT  IS  'LETTER 
(t 

(putprop  pixname  'flat  'sec  class) 

(putprop  pixname  flat  char  'letter)))) 

ELSE.  IF  THE  SLOPE  IS  LESS  THAN  THE  UNSAFE  SLOPE  LIMIT  THEN 
THE  PRIMARY  CLASSIFICATION  IS  SAFE  SLOPE 
-AND  SAFE  SLOPE  IS  THE  LETTER 
((<  slope  unsafe  slope  limit) 

(putprop  pixname  safe  slope  char  'letter) 

(putprop  pixname  'safe  slope  ’priclass)) 

ELSE  EVERYTHING  ELSE  IS  UNSAFE  SLOPE  PRIMARY  CLASSIFICATION 
;-AND  UNSAFE  SLOPE  IS  THE  LETTER' 

(t 

(putprop  pixname  'unsafe  slope  'pri  class) 
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(putprop  pixname  unsafe  slope  char  letter)) 

:AND  MAX  AND  MIN  CURVATURE  IS  NEGATIVE 
.-PEAK  IS  SECONDARY  CLASSIFICATION 

(cond  ((and  (<  El  (minus  curv  threshold)) 

(<  E2  (minus  curv  threshold))) 

(putprop  pixname  ’peak  ’secclass)) 

:0R  MAX  AND  MIN  CURVATURE  IS  POSITIVE 
;-PIT  IS  SECONDARY  CLASSIFICATION 
((and  (>  El  curv_threshold) 

(>  E2  curv  threshold)) 

(putprop  pixname  ’pit  ’sec_dass)) 

:OR  MAX  CURVATURE  IS  NEGATIVE  WHILE  MIN  CURVATURE  IS  APPROX  0 
;-RIDGE  IS  SECONDARY  CLASSIFICATION 
((and  (<  El  (minus  curv  threshold)) 

(<  (abs  E2)  curv  threshold)) 

(putprop  pixname  ’ridge  ’secclass)) 

:OR  MAX  CURVATURE  IS  POSITIVE  WHILE  MIN  CURVATURE  IS  APPROX  0 
:-RAVINE  IS  SECONDARY  CLASSIFICATION 
((and  (>  El  curv  threshold) 

(<  (abs  E2)  curv_threshold)) 

(putprop  pixname  ’ravine  ’secclass)) 

;OR  IF  THE  MAX  CURVATURE  IS  NEGATIVE  WHILE  THE  MIN  CURVATURE 
:POSITIVE 

;-SADDLE  IS  SECONDARY  CLASSIFICATION 
((and  (<  El  (minus  curv  threshold)) 

(>  E2  curv  threshold)) 

(putprop  pixname  ’saddle  ’sec  class)) 

OR  IF  THE  MAX  CURVATURE  IS  POSITIVE  WHILE  THE  MIN  CURVATURE 
;IS  NEGATIVE 

;-PASS  IS  SECONDARY  CLASSIFICATION 
((and  (>  El  curv  threshold) 

(<  E2  (minus  curv  threshold))) 

(putprop  pixname  ’pass  ’sec  class)) 

:  EVERYTHING  ELSE  IS 

;-FLAT  TERRAIN  IS  SECONDARY  CLASSIFICATION 
(‘ 

(putprop  pixname  ’flat  ’sec  class)))))) 

**  "draw  rclass  map"  is  the  function  to  call  to  create 
**  a  real  classification  map,  draw  the  class  map  base 
**  to  a  bit  map.  and  then  draw  the  classification 
**  map  to  the  bit  map. 

(defun  draw  rclass  map  (bmd  fd) 

(prog  (pi  pj  map  x  cmx  cmy) 

(make  rclass  map) 

( draw  class  map  base  bmd) 

(setq  pi  0) 

(setq  pj  0) 

(setq  map  rclass_map) 
loop 

(setq  x  (car  map)) 
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(setq  map  (cdr  map)) 

(setq  cmx  (—  xOcminit  (*  14  pi))) 

(setq  cmy  (—  yOcminit  (*  14  pj ) ) ) 

(cond  ((=  (get  x  ’letter)  7) 

(paint  safe  slope  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’U) 

( paint  unsafe  slope  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’E) 

(paint  edge  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’-) 

(paintflat  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’P) 

(paint  peak  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’L) 

(paint  pit  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’*) 

(paint  ridge  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’v) 

(paintravine  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’Y) 

(paint  pass  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’S) 

(paint  saddle  bmd  cmx  cmy))) 

(cond  ((and  (  =  pi  numacross) 

(=  pj  numdown)) 

(bmdisp  fd  COP YR ASTER  bmd  0  0  0  0  wwidth  w  height  0) 
(bmadres  bmd  VT  ABSOLUTE) 

(return)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1+  pj)) 

(go  loop)) 

(' 

(setq  pi  (1-e  pi)) 

(go  loop))))) 

**  "draw  nclass  map"  is  the  function  to  call 
**  to  create  a  noisy  classification  map, 

**  draw  the  class  map  base  to  a  bit  map,  and  then 
**  draw  the  classification  map  to  the  bit  map. 

(defun  draw  nclass  map  (bmd  fd) 

(prog  (pi  pj  map  x  cmx  cmy) 

(makenclassmap) 

(draw _class_map_base  bmd) 

(setq  pi  0) 

(setq  pj  0) 

(setq  map  nclass  map) 
loop 

(setq  x  (car  map)) 

(setq  map  (cdr  map)) 

(setq  cmx  (+  xOcminit  (*  14  pi))) 
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(setq  cmy  (—  yOcminit  (’  14  pj ) ) ) 

(cond  ((=  (get  x  ’letter)  ’;) 

(paint  safe  slope  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’U) 

(paint  unsafe  slope  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’E) 

(paintedge  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’-) 

(paintflat  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’P) 

(paint  peak  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’L) 

(paintpit  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’  *) 

(paint  ridge  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’v) 

(paint  ravine  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’Y) 

(paint  pass  bmd  cmx  cmy)) 

((=  (get  x  ’letter)  ’S) 

(paint  saddle  bmd  cmx  cmy))) 

(cond  ((and  (  =  pi  numacross) 

(=  pj  numdown)) 

(bmdisp  fd  COPYR ASTER  bmd  0  0  0  0  wwidth  wheight  0) 
(bmadres  bmd  VT  ABSOLUTE) 

(return)) 

((=  pi  numacross) 

(setq  pi  0) 

(setq  pj  (1+  pj)) 

(go  loop)) 

(t 

(setq  pi  (1+  pi)) 

(go  loop))))) 

**  "drawl"  draws  the  first  display,  that  being 
**  the  real  terrain 

(defun  drawl  () 

(initialize) 

(make  real _terrai n ) 

(draw  terrain  ’real  terram  (get(car  windows)’bmd) 

(get(car  windowsj’fdnum)) 

(copy  bm  tofile  (get  (car  windows)  ’bmd)) 

( make  noisy  terrain ) 

(draw  terrain  ’noisy  terrain  (  get  (cadr  windows)  ’bmd) 

(  get  (cadr  windows)  ’fdnum)) 

(copy  bm  tofile  (get  (cadr  windows)  ’bmd)) 

(draw  rclass  map  (get  (caddr  windows)  ’bmd) 

(get  (caddr  windows)  ’fdnum)) 

(copy  bm  tofile  (get  (caddr  windows)  ’bmd)) 

(draw _nc lass  map  (get  (cadddr  windows)  ’bmd) 

(get  (cadddr  windows)  ’fdnum)) 
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(copy_bm_tofile  (get  (cadddr  windows)  ’bmd))) 

**  "redraw"  redraws  a  new  real  terrain  into  window  1 

(defun  redraw  () 

(prog  (bm  fd) 

(setq  bm  (get  (car  windows)  ’bmd)) 

(setq  fd  (get  (car  windows)  ’fdnum)) 

(bmciear  bm  0  0  wwidth  wheight  VTWHITE) 

( make  real  terrain ) 

(draw  terrain  ’realterrain  bm  fd) 

(copybmtofile  bm) 

(setq  bm  (get  (cadr  windows)  ’bmd)) 

(setq  fd  (get  (cadr  windows)  ’fdnum)) 

(bmciear  bm  0  0  wwidth  wheight  VT_WHITE) 

( make  noisy  terrain ) 

(draw  terrain  ’noisy  terrain  bm  fd) 

(copy  bm  tofile  bm) 

(setq  bm  (get  (caddr  windows)  ’bmd)) 

(setq  fd  (get  (caddr  windows)  ’fdnum)) 

(bmciear  bm  0  0  wwidth  wheight  VT  WHITE) 
(drawrclassmap  bm  fd) 

(copy  bm  tofile  bm) 

(setq  bm  (get  (cadddr  windows)  ’bmd)) 

(setq  fd  (get  (cadddr  windows)  ’fdnum)) 

(bmciear  bm  0  0  wwidth  wheight  VT  WHITE) 

(draw  nclass  map  bm  fd) 

(copy  bm  tofile  bm))) 

1 

;  **  "keep  running"  starts  the  program 

**  with  a  call  to  "drawl"  then  it  checks 
;  **  to  see  if  it  has  completed  the  required 

**  number  of  runs.  If  it  has  it  closes  all 
**  the  windows  and  exits  lisp.  If  not  it 
**  keeps  running. 

(defun  keep  running  () 

(prog  () 

(drawl) 
loop  10 

(cond  ((=  fiienum  (-i-  (*  4  numruns)  filebase)) 

(finish) 

(exit)) 

(‘ 

(redraw) 

(go  looplO))))) 

;  **  "finish"  closes  all  the  open  windows  and 

;  **  deallocates  the  bitmap  memory  space. 

1 

(defun  finish  () 

(prog  (windowlist  x) 
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(setq  windowlist  windows) 
loop 

(setq  x  (car  windowlist)) 

(setq  windowlist  (cdr  windowlist)) 

(clw  (get  x  ’fdnum)) 

(bmdeal  (get  x  ’bmd)) 

(cond  ((  not  (null  windowlist))  (go  loop))))) 

******************************************************** 
*  * 

|  *  THE  MAIN  PROGRAM  * 

*  * 

******************************************************** 

;  **  Define  window  one 

(def  window  ’wl  wlx  wly 

wwidth  wheight  "TERRAIN") 

**  Define  window  two 

(defwindow  ’w2  wlx  (plus  wly  ysep  wheight) 
wwidth  wheight  "NOISY  TERRAIN") 

**  Define  window  three 

(def  window  'w3  (plus  wlx  xsep  wwidth)  wly 
wwidth  wheight  "CLASSIFICATION  MAP") 

;  **  Define  window  four 

(def  window  ’w4  (plus  wlx  xsep  wwidth)  (plus  wly  ysep  wheight) 
wwidth  wheight  "NOISY  CLASSIFICATION  MAP") 

**  Start  the  simulation 

(keep_running) 

it************************************************* 

*  * 

;  *  THE  END 

*  * 

*************************************************** 
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